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(ABSTRACT) 


Highly  agile,  hover  capable  flapping  wing  flight  is  a  relatively  new  area  of  study  in  engi¬ 
neering.  Researchers  are  looking  to  flapping  flight  as  a  potential  source  for  the  next  gener¬ 
ation  of  reconnaissance  and  surveillance  vehicles.  These  systems  involve  highly  complicated 
physics  surrounding  the  flapping  wing  motion  and  unusual  characteristics  due  to  a  hover  re¬ 
quirement  not  normally  associated  with  conventional  aircraft.  To  that  end  this  study  focuses 
on  examining  the  various  models  and  physical  parameters  that  are  considered  in  various 
other  studies.  The  importance  of  these  models  is  considered  through  their  effect  on  the  trim 
and  stability  of  the  overall  system.  The  equations  of  motion  are  modeled  through  a  quasi  co¬ 
ordinate  Lagrangian  scheme  while  the  aerodynamic  forces  are  calculated  using  quasi-steady 
potential  flow  aerodynamics.  Trim  solutions  are  calculated  using  periodic  shooting  for  sev¬ 
eral  different  conditions  including  hover,  climb,  and  forward  flight.  The  stability  of  the  trim 
is  calculated  and  examined  using  stroke-averaged  and  Floquet  theory.  Inflow  and  viscous  ef¬ 
fects  are  added  and  their  effects  on  trim  and  stability  examined.  The  effects  of  varying  hinge 
location  and  the  inclusion  of  stroke  deviation  in  the  wing  kinematics  are  also  explored.  The 
stroke-averaged  system  was  not  found  to  be  a  direct  replacement  for  the  periodic  system  as 
the  stability  was  different  for  the  two  systems.  Inflow  and  viscosity  were  found  to  have  large 
effects  on  the  stability  of  the  system  and  models  accounting  for  the  two  should  be  included 
in  future  flight  dynamic  models. 
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Chapter  1 
Introduction 


1.1  Motivation 

The  Defense  Advanced  Research  Project  Agency  (DARPA)  defines  a  Micro  Air  Vehicle 
(MAV)  as  an  aircraft  limited  in  all  of  its  dimensions  to  15  cm  [4].  This  size  limitation 
is  meant  to  allow  the  vehicle  to  operate  in  a  confined  space  such  as  those  presented  by  an 
urban  canyon  or  building  interior.  While  any  type  of  aircraft  could  theoretically  be  adap- 
ated  to  fit  the  MAV  definition,  rotor  craft  and  flapping  wing  designs  are  particularly  well 
adapted  for  the  high  agility  and  maneverability  requirements  of  this  difficult  operating  envi¬ 
ronment.  Further  more,  for  concealment  in  plain  sight,  flapping  wing  designs  are  particularly 
well-suited  in  that  they  can  mimic  designs  we  see  in  nature  and  thereby  have  the  potential 
to  hide  in  plain  sight.  Of  particular  interest  in  demonstrating  these  capabilities  are  recent 
developments  in  highly  agile,  hover  capable,  flapping-wing  aircraft  that  mimic  small  fliers  in 
nature  such  as  Aerovironment’s  Nano  Hummingbird  [5]. 

This  particular  study  focuses  on  the  difficulty  of  determining  what  models  are  most  impor¬ 
tant  to  consider  when  trying  to  accurately  represent  flapping  wing  dynamics.  Consideration 
is  given  to  eventual  control  synthesis  and  calcuation  of  dynamic  metrics  such  as  performance, 
stability,  maneuverability,  and  gust  response.  To  this  end,  various  models  such  as  viscous 
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forces  and  inflow,  and  physical  parameters,  such  as  deviation  and  hinge  location  are  con¬ 
sidered.  These  models  are  incorporated  into  an  otherwise  quasi-steady  aerodynamic  model 
(integrated  over  the  span  using  blade  element  theory)  which  is  coupled  with  a  rigid-body 
quasi-coordinate  Lagrangian  scheme. 


1.2  Literature  Review 

The  literature  review  consists  of  three  sections.  The  first  deals  with  flight  dynamics  and 
equations  of  motion  with  an  emphasis  on  those  methods  that  have  been  applied  to  flapping 
wing  systems.  The  second  section  deals  with  the  various  aerodynamic  models  that  have  been 
used  on  flapping  wings.  The  third  deals  with  stability  and  control. 

1.2.1  Aerodynamics  of  Flapping  Wing  Systems 

Fry  et  al.  [6]  examined  the  kinematics  and  aerodynamics  of  Drosophila  flight,  specifically 
rapid  maneuvers,  using  high-speed  infrared  photography.  The  captured  kinematics  were 
played  through  a  dynamically  scaled  robot  in  order  to  measure  the  produced  aerodynamic 
forces.  The  maneuvers  in  question  are  described  as  saccades,  fast  turns  that  can  involve 
as  much  as  ninety  degree  turns  in  less  than  fifty  micro  seconds.  Results  showed  that  the 
mean  forces  generated  by  the  fruit  fly  were  fairly  constant  with  respect  to  the  body  axes. 
Instead,  during  maneuvers,  the  body  orientation  was  changed  in  order  to  facilitate  overall 
acceleration  (along  with  very  slight  changes  in  wing  motion).  The  paper  also  examines  the 
assumption  that  smaller  animals  must  overcome  primarily  viscous  friction  forces  while  larger 
animals  are  dominated  by  inertial  forces  when  accelerating  or  turning.  Through  analysis  of 
the  time  course  data  for  the  torque  produced  by  the  insect,  it  was  concluded  that  the  very 
small  fruit  fly  was  not  dominated  by  viscous  friction  forces  but  instead  by  inertial  forces.  It 
is  also  stated  that  while  the  research  was  conducted  on  the  fruit  fly,  the  results  are  relevant 
for  nearly  all  insects  in  that  inertial  forces  only  become  more  important  as  size  increases  and 
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the  turns  observed  were  banked  turns  similar  to  those  observed  in  larger  insect  species. 

Continuing  this  work,  in  Fry  [7]  the  body  and  wing  kinematics  of  Drosophila  flies  were 
characterized  through  the  use  of  high-speed  video,  lasers,  and  a  dynamically  scaled  robot. 
The  motion  was  captured  using  high-speed  cameras  focused  on  a  0.5  cm 3  area  in  a  control 
volume  lit  by  infrared  LEDs  (which  are  not  visible  to  the  flies  and  so  do  not  disrupt  their  be¬ 
havior).  The  robot  was  then  used  to  ‘replay’  motions  captured  by  the  video  and  measure  the 
aerodynamic  forces  caused  by  these  motions.  A  U-shaped  motion  was  observed  in  the  wings 
of  the  hovering  flies  with  no  observed  ‘clap  and  fling’  motion.  Quasi-steady  aerodynamics 
were  found  to  be  able  to  account  for  nearly  all  of  the  behavior  involved  in  hovering  except 
that  there  are  discrepancies  between  measured  and  model  data  that  indicate  the  presence 
of  unsteady  mechanisms.  The  paper  also  details  the  comparison  between  tethered  and  un¬ 
tethered  flight  of  the  flies  and  how  the  stroke  pattern  is  severely  disrupted  by  the  tethering. 
This  refutes  the  historical  methods  of  using  tethered  flight  as  a  proxy  for  hovering  flight, 
instead  showing  that  the  tethering  makes  for  different  stroke  patterns  and  a  significant  pitch 
down  moment.  Fry  proposes  that  this  is  an  effect  of  sensory  disruption.  They  conclude  by 
describing  the  balance  necessary  between  flight  forces  for  stable  hovering  and  the  challenge 
this  produces  with  two  high  frequency  wings.  They  propose  that  the  complex  nature  of  the 
still  poorly  understood  hinge  may  hold  the  key  to  permitting  a  limited  number  of  muscles  to 
make  the  subtle  changes  in  strokes  necessary  to  control  insect  flight.  A  table  providing  the 
parameters  and  power  numbers  measured  in  the  experiment  is  also  provided  as  an  important 
resource. 

Singh  and  Chopra  [8]  examine  the  aerodynamics  of  insect-based,  biomimetic,  flapping 
wings  in  hover.  A  biomimetic  flapping  wing  test  mechanism  was  used  to  imitate  insect 
wing  flapping.  The  mechanism  was  a  passive-pitch  bi-stable  mechanism  controlled  by  a 
speed  controller,  pulse  generator,  and  brushless  motor.  Due  to  vibration  and  inertial  loads, 
mounting  the  apparatus  on  a  load  cell  did  not  provide  very  meaningful  results.  Instead,  load 
cells  made  of  strain  gages  were  mounted  at  the  base  of  the  wings  to  measure  the  moments 
caused  by  the  wings.  The  motion  of  the  wings  (pitch  angle)  was  measured  using  a  Hall  sensor 
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(a  magnet  moving  past  the  Hall  effect  sensor  causes  a  relative  charge).  Inertial  loads  were 
measured  by  running  the  mechanism  in  a  vacuum  chamber  (90%  vacuum).  Low  frequency 
tests  were  run  for  wings  I-III  and  results  reported  for  wings  II  and  III.  Wing  I  and  II  were  of 
the  same  form  (wing  II  had  thicker  Mylar)  which  was  an  aluminum  frame  around  all  edges 
of  the  frame  and  the  wing  surface  itself  made  of  Mylar.  Wing  III  was  similar  but  the  frame 
did  not  continue  around  the  trailing  edge,  but  instead  only  around  the  leading  edge  and  the 
root.  High  frequency  results  were  reported  for  wing  III,  and  several  tests  were  also  carried 
out  for  a  series  of  lightweight  wings  that  used  carbon  fiber  frames  rather  than  aluminum 
frames.  Results  were  also  given  for  several  of  these  light  wings.  The  analytical  model  for  the 
forces  accounted  for  translational  and  rotational  circulation,  the  leading-edge  vortex  effects, 
noncirculatory  forces  based  on  thin  airfoil  theory,  starting  vortex  effects,  and  shed  wake  and 
tip  vortex  effects.  At  low  frequency  wing  III  showed  the  most  thrust  (pitched  about  its  20% 
chord  location)  and  showed  more  thrust  at  a  pitch  angle  of  45  degrees  rather  than  30  degrees. 
When  the  frequency  was  increased  to  try  and  increase  thrust  the  mass  of  the  wing  was  shown 
to  have  a  significant  effect  on  the  maximum  possible  frequency  to  which  the  mechanism  could 
drive  the  wing.  It  was  also  noted  that  the  thrust  showed  a  sudden  drop  at  a  frequency  of 
11.6  Hz  (wing  III).  All  wings  showed  an  eventual  drop  at  high  frequency,  possibly  due  to  a 
loss  in  noncirculatory  thrust  caused  by  elasticity.  It  was  also  noted  that  using  a  softer  spring 
to  control  the  passive  pitching  caused  a  larger  pitch  variation  that  resulted  in  a  larger  thrust 
at  a  smaller  power  consumption.  Singh  and  Chopra  propose  that  this  passive  mechanism 
merits  further  study  due  to  its  mechanical  simplicity  and  ease  of  implementation. 

Dickinson  et  al.  [2]  discuss  the  interaction  of  three  interactive  mechanisms  that  enhance 
aerodynamic  performance.  The  three  mechanisms  are  delayed  stall,  rotational  circulation, 
and  wake  capture.  Delayed  stall  occurs  during  the  translational  portions  of  the  stroke  while 
rotational  circulation  and  wake  capture  occur  during  the  stroke  reversal.  These  rotational 
mechanisms  can  modulate  the  direction  and  magnitude  of  flight  forces  during  steering  ma¬ 
neuvers.  When  insect  wings  are  placed  in  a  wind  tunnel  and  tested  over  the  range  of  ve¬ 
locities  normally  encountered  the  measured  forces  are  much  lower  than  those  required  for 
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flight.  They  also  discuss  one  of  the  first  unsteady  effects  to  be  identified:  clap  and  fling. 
This  phenomenon  precedes  pronation  and  hastens  the  development  of  circulation  during  the 
downstroke.  Clap  and  fling  is  not  present  in  all  insects  and  thus  cannot  be  the  solution  to 
the  lower  than  required  forces  presented  in  the  wind  tunnel  tests.  Delayed  stall  represents 
a  flow  structure  on  the  leading  edge  of  a  wing  that  can  generate  circulatory  forces  greater 
than  those  experienced  under  steady-state  conditions  at  high  angles  of  attack.  The  3D  force 
coefficients  are  greater  than  those  of  the  transient  2D  steady-state  values  which  confirm  the 
presence  of  some  transient  unsteady  effect  such  as  delayed  stall.  Although  the  translational 
values  with  the  delayed  stall  effects  match  the  magnitude  of  the  measured  force  near  the 
middle  of  each  half  stroke  they  do  not  match  the  values  during  force  reversal.  A  possible 
explanation  for  the  peaks  at  the  end  of  each  half  stroke  could  be  a  combination  of  rotational 
effects  and  wake  capture.  Rotational  effects  are  similar  to  those  effects  experienced  by  a  spin¬ 
ning  baseball;  air  is  pulled  within  the  boundary  layer  and  serves  as  a  source  of  circulation. 
The  wing  must  rotate  to  be  in  the  proper  position  for  each  stroke.  This  rotation,  however, 
can  be  altered  to  change  the  reaction  forces  caused  by  the  rotational  effects.  If  the  wing  flips 
early  (before  changing  direction)  then  the  resulting  force  should  be  upward.  If  the  flip  is  late 
then  the  resulting  force  would  be  a  downward  force.  If  the  flip  spans  both  half  strokes  then 
the  resulting  downward  and  upward  forces  cancel.  Therefore  proper  timing  of  wing  rotation 
can  add  extra  lift  beyond  that  of  just  delayed  stall.  Delayed  stall,  however,  does  not  account 
for  the  large  positive  transient  that  develops  immediately  after  the  wing  changes  direction. 
The  researchers  propose  that  wake  capture  is  the  other  rotational  effect  that  makes  up  the 
difference  between  the  translational  effects  and  the  rotational  effects.  In  wake  capture  the 
wing  benefits  from  the  shed  vorticity  of  the  previous  stroke.  The  vorticity  generated  by 
the  previous  stroke  increases  the  effective  velocity  of  the  next  stroke  and  increases  the  force 
production.  To  support  this  hypothesis,  the  authors  tested  the  force  produced  by  halting  the 
wing  at  the  end  of  the  half  stroke.  The  results  showed  that  wing  generates  force  for  several 
hundred  milliseconds  following  the  end  of  the  half  stroke.  Again,  the  rotation  timing  has  a 
similar  effect  to  that  of  delayed  stall.  The  authors  conclude  that  the  wing  should  rotate  in 
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advance  of  the  reversal  so  that  it  generates  circulatory  forces  at  the  end  of  each  stroke,  and 
increases  the  strength  of  the  wake  and  ensures  proper  orientation  to  take  advantage  of  the 
wake.  Insects  alter  this  timing  in  order  to  aid  in  control.  For  instance  the  rotation  timing 
on  one  wing  can  be  altered  to  create  more  drag  and  aid  in  turning. 

Zbikowski  [9]  discusses  advances  in  fluid  dynamics  as  applied  to  insects  and  micro  air 
vehicles  (MAVs).  These  new  concepts  are  applied  to  two  methods  of  aerodynamic  modeling 
of  an  insect  like  wing  in  hover.  One  of  the  unsteady  aerodynamic  phenomena  Zbikowski 
discusses  is  a  spiraling  leading-edge  vortex.  This  vortex  is  bound,  remaining  on  the  wing 
during  the  half  cycle  throughout  pitching,  plunging,  and  sweeping  before  being  shed  at  the 
end  of  the  stroke  when  the  wing  flips.  This  is  very  similar  to  Dickinson’s  delayed  stall 
effect.  The  other  effects  included  those  of  wing  pitching,  plunging,  and  sweeping  present  at 
all  times  and  wing  interactions  with  its  own  wake  (Dickinson’s  wake  capture  effect).  The 
flow  was  assumed  incompressible,  low  Reynolds  number,  and  laminar  with  a  rigid,  thin  wing 
of  symmetric  section.  Zbikowski  argues  that  these  assumptions  are  well  supported  except 
for  lack  of  flexibility  in  the  wing  which  he  adopts  for  simplicity.  He  proposes  that  insects 
deliberately  force  separation  at  the  leading  edge  to  gain  vortical  lift  and  shedding  the  vortex 
at  the  end  of  the  half  stroke  using  the  sudden  wing  flip.  The  pitching,  plunging,  and  sweeping 
nature  of  the  wing  are  similar  to  aspects  of  helicopter  blades  but  with  the  critical  difference 
that  MAV  wings  consistently  flap  with  angles  of  attack  greater  than  20  degrees  without  having 
the  vortex  shed  and  losing  lift.  In  insect  like  flapping  the  separation  occurs  at  the  beginning  of 
the  motion  and  generates  stable  lift  (from  the  vortex)  throughout  the  half  cycle.  The  essence 
of  Zbikowski’s  framework  is  to  account  for  both  the  leading-edge  vortex  and  the  other  part  of 
the  flow.  The  non-vortical  part  of  the  flow  was  treated  with  unsteady  aerodynamics  similar 
to  those  used  in  helicopter  blades  with  attached  flow,  while  the  bound  leading-edge  vortex 
was  considered  a  non-linear  phenomenon  and  had  to  be  treated  with  methods  different  from 
the  thin  airfoil  theory.  Using  the  circulation  approach  allows  for  the  use  of  just  a  single 
framework  since  both  the  linear  and  nonlinear  contributions  can  be  handled  consistently. 
Using  the  velocity  potential  necessitates  using  a  second  correction  to  deal  with  the  nonlinear 
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components  (in  this  case  Polhamus’s  leading-edge  suction  analogy).  The  circulation  approach 
uses  a  methodology  based  on  linear,  thin-airfoil  theory  applied  to  unsteady  wing  motions. 
The  linear  portions  of  the  aerodynamics  are  handled  with  classical  von  Karman-Sears  Theory 
and  the  nonlinear  portions  with  McCune-Tavares  nonlinear  extensions.  The  velocity  potential 
approach  is  based  upon  the  theories  of  Theodorsen,  Wagner,  and  Garrick  supplemented  by 
the  Polhamus  leading  edge  suction  analogy.  Zbikowski  argues  that  the  simplifications  made 
in  his  paper  are  necessary  now  in  order  to  further  engineering  design  so  that  fundamental 
analyses  will  be  furthered,  leading  to  production  of  a  better  design. 

Nagai  et  al.  [10]  detail  experimental  and  numerical  studies  on  the  aerodynamic  character¬ 
istics  of  an  insect  in  forward  flight  (specifically  the  flapping  wing).  Much  like  in  Dickinson’s 
work,  the  unsteady  aerodynamic  forces  and  flow  patterns  were  measured  using  a  scaled  me¬ 
chanical  model  but  in  this  case  in  a  water  tunnel  instead  of  mineral  oil.  The  wing  and 
motions  were  based  on  the  flapping  wing  of  a  bumblebee  and  the  results  were  compared  to 
a  3D  Navier-Stokes  code.  The  aerodynamic  mechanisms  of  the  wing,  such  as  delayed  stall, 
rotational  effect,  and  wake  capture,  were  examined  in  detail  for  forward  flight.  The  motions 
played  through  the  model  were  ‘simple  trapezoidal-type’  motions  based  on  data  from  Dudley 
and  Ellington  [11],  The  mechanism  was  similar  to  that  of  Dickinson’s  in  that  the  mechanism 
provided  the  three  degrees  of  freedom  and  in  the  cases  of  forward  flight  the  water  was  moved 
through  the  tunnel  at  0.02  to  0.2  m/s.  Because  a  bumblebee  moves  its  hind  and  forewings 
together  and  not  independently,  the  planform  chosen  for  this  wing  encompassed  both  and 
was  made  of  acrylic  surrounded  by  an  aluminum  frame.  Forces  were  measured  at  the  root  by 
strain  gages  in  the  typical  fashion.  For  flow  visualization  the  wing  was  slightly  different  in 
that  it  had  no  aluminum  frame  so  that  light  could  pass  through  the  entire  wing  and  it  was 
made  thicker  to  retain  resistance  to  bending.  Comparisons  between  the  two  wings  showed  2% 
or  less  difference  in  aerodynamic  characteristics.  A  laser  shining  through  the  side  wall  of  the 
water  tunnel  was  used  as  a  light  source  and  two  cameras  placed  below  the  tunnel  recorded 
the  flow  patterns.  The  three  dimensional  Navier-Stokes  code  has  been  validated  by  Isogai 
et  al.  [12,  13]  in  relation  to  the  time  histories  (hovering  flight)  of  aerodynamic  forces  for  a 
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dragonfly.  Three  cases  were  chosen  for  validation  in  this  experiment,  a  course  grid,  a  fine 
grid,  and  a  fine  time-step  case.  All  three  cases  were  in  good  agreement  with  the  experimental 
results.  Flow  patterns  were  also  compared  between  the  experimental  and  numerical  results 
and  again  agreement  between  the  two  was  good.  There  were  small  quantitative  differences 
throughout  the  comparisons  that  were  attributed  to  the  differences  in  shape  between  the 
numerical  and  experimental  wings.  The  experimental  wing  had  a  thickness  as  described  pre¬ 
viously  with  a  rectangular  leading  edge  while  the  numerical  model  had  a  sharp  leading  edge 
with  zero  thickness.  This  leads  to  a  larger  leading  edge  vortex  in  the  experimental  case  than 
in  the  numerical  one  which  may  account  for  the  difference  in  quantitative  results. 

Sane  and  Dickinson  [14]  compared  instantaneous  force  measurements  from  a  dynamically- 
scaled  mechanical  flapping  wing  to  quasi-steady  estimates  in  order  to  study  rotational  forces 
for  a  range  of  angular  velocities  and  axes  of  rotation.  It  was  found  that  the  rotational 
coefficient  varied  as  expected  with  axis  of  rotation  but  also  varied  with  angular  velocity  (not 
expected  based  on  theoretical  models).  The  standard  quasi-steady  model  was  then  modified 
to  include  rotational,  translational,  and  added-mass  effects.  By  subtracting  this  new  quasi¬ 
steady  estimate  from  the  measured  forces  the  aerodynamic  effects  due  to  wake  capture  were 
isolated.  This  paper  also  gives  a  brief  description  of  what  makes  up  each  of  the  effects  added 
to  the  standard  quasi-steady  estimate. 

Peters  et  al.  [15]  extend  classical  thin  airfoil  theory  by  allowing  for  large  frame  motion 
of  the  wing  and  small  deformations  (expanded  into  Chebychev  polynomials)  relative  to  large 
frame  motion  (in  particular  for  rotorcraft  analysis).  The  theory  itself  is  formulated  in  terms 
of  the  deflections  and  generalized  forces  within  the  frame  and  is  based  on  transformations  of 
the  traditional  potential  flow  equations.  The  equations  are  broken  down  into  mass,  damping 
and  stiffness  matrices,  making  coupling  with  finite  elements  or  modal  analysis  convenient. 
Also  the  unsteady  components  of  the  aerodynamics  are  encompassed  in  a  completely  separate 
wake  model  making  the  theory  flexible  in  that  it  can  be  changed  to  match  the  physics  rather 
than  relying  on  a  constant  unsteady  model  such  as  the  flat  wake  found  in  Theodorsen’s 
theory.  This  allows  for  a  theory  somewhere  in  between  quasi-steady  blade  element  theories 
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and  CFD  analyses  in  terms  of  accuracy  and  computation  time.  Peters  et  al.  validated  the 
theory  against  those  of  Theodorsen  [16],  Garrick  [17],  Loewy  [18],  Issacs  [19,  20],  Greenberg 
[21],  Wagner  [22],  and  von  Karman  [23]  as  well  as  experimental  cases. 

Pesavento  and  Wang  [24]  solved  the  two  dimensional  Navier-Stokes  equations  in  relation 
to  a  piece  of  free  falling  paper  at  Reynolds  number  near  103.  In  contrast  to  an  airfoil 
governed  by  the  Kutta  Joukowsld  condition  which  normally  experiences  lift  proportional  to 
velocity  squared,  the  lift  was  found  to  be  dominated  by  the  product  of  the  linear  and  angular 
velocities.  In  comparing  the  two  models,  the  Joukowsld  model  was  modified  to  incorporate 
a  lift-drag  polar  but  this  still  did  not  exhibit  the  cusplike  turning  points  with  elevation  of 
the  center  of  mass.  This  difference  was  further  investigated  through  decomposition  of  forces 
into  added-mass  and  circulation  terms  and  the  Navier-Stokes  equations  are  found  to  provide 
an  extra  instantaneous  forcing  term  in  the  circulation  due  to  the  angular  velocity. 

Anderson  et  al.  [25]  continued  the  earlier  work  in  [24],  this  time  experimentally  quantify¬ 
ing  the  trajectories  of  falling  plates.  The  trajectories  were  captured  using  high-speed  digital 
video  in  order  to  determine  instantaneous  accelerations  from  which  to  deduce  the  fluid  forces. 
These  values  were  compared  to  direct  solutions  of  2D  Navier-Stokes  equations.  Again  the 
fluid  circulation  was  found  to  be  dominated  by  a  rotational  term  proportional  to  the  angular 
velocity  of  the  plate,  unlike  the  lift  typically  found  for  an  airfoil  at  an  angle  of  attack.  The 
torque  was  also  found  to  be  small  in  comparison  to  that  normally  found  on  an  airfoil  at  a 
fixed  angle  of  attack.  By  varying  Reynolds  number,  dimensionless  moment  of  inertia,  and 
the  thickness  to  width  ratio  of  the  plate,  the  dynamics  varied  from  fluttering,  tumbling,  and 
chaotic  motion.  Variations  in  initial  conditions  provided  a  period  of  brief  transients  followed 
by  periodic  fluttering.  These  results  matched  well  with  the  numerical  results  previously 
described  in  [24], 

Berman  and  Wang  [26]continued  previous  work  done  by  Wang  and  others  in  [24,  25]  by  in¬ 
vestigating  aspects  of  hovering  insect  flight.  Optimal  wing  kinematics  that  minimized  power 
usage  while  constrained  to  provide  enough  lift  for  hover  (maintain  time-averaged  constant 
altitude  over  a  flapping  period)  were  determined  using  a  hybrid  optimization  that  combines 
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aspects  of  genetic  and  gradient  based  optimizers.  The  wings  examined  in  this  study  were 
considered  as  rigid  bodies  with  three  rotational  degrees  of  freedom  and  are  based  on  the 
wings  of  a  fruit  fly,  bumblebee,  and  hawkmoth.  In  order  to  allow  for  the  many  cost  func¬ 
tion  evaluations  and  detailed  sensitivity  analysis  that  go  into  the  optimization  procedure  the 
quasi  steady  aerodynamics  described  in  [24]  were  utilized.  The  results  of  the  analysis  pro¬ 
vided  kinematics  which  were  qualitatively  and  quantitatively  similar  to  observed  data  ([7,  27] 
among  others).  The  data  also  showed  a  constant  leading  edge  throughout  the  stroke  (trailing 
edge  never  lead)  being  advantageous  due  to  interactions  between  inertial  and  aerodynamic 
power. 

Stanford  et  al.  [3]  examined  the  effects  of  several  parameters  on  aeroelastic  hovering 
motions  of  a  flexible  flapping  wing.  The  parameters  used  in  this  study  were  wing  shape, 
structural  composition,  and  the  kinematic  motions,  while  the  criteria  for  performance  was 
the  peak  power  required  for  the  stroke  subject  to  kinematic  constraints.  Because  this  work 
was  completed  at  the  partner  organization  Air  Force  Research  Laboratory  (AFRL)  the  aero¬ 
dynamics  of  this  study  were  of  particular  interest  for  comparison  and  validation.  The  aerody¬ 
namics  consisted  of  a  quasi-steady  blade  element  model  and  a  higher-order  Navier  Stokes  flow 
solver  (OVERFLOW-2.1).  The  quasi  steady  aerodynamics  were  based  on  those  of  Berman 
and  Wang  [26]  and  differ  from  those  in  this  study  due  to  the  inclusion  of  viscous  forcing 
terms  as  well  as  empirically  determined  forcing  coefficients.  With  these  terms  included  the 
quasi-steady  aerodynamics  were  able  to  show  good  agreement  with  the  CFD  model  despite 
the  large  difference  in  computational  time  and  model  fidelity. 

Oppenheimer  et  al.  [28]  developed  an  aerodynamic  model  from  blade  element  theory 
(not  including  unsteady  effects)  based  on  the  work  of  Sane  and  Dickinson  [14],  The  vehicle 
considered  in  the  work  is  similar  to  the  Harvard  RoboFly  except  that  the  wings  are  indepen¬ 
dently  actuated  using  two  different  piezo  electric  actuators.  The  blade  element  theory  is  used 
to  calculate  instantaneous  and  stroke-averaged  forcing  and  moments  for  use  in  control  of  the 
aircraft.  The  position  of  the  wing  is  controlled  using  oscillators  that  change  once  per  wing 
beat  cycle  so  as  not  to  invalidate  the  stroke-averaged  model.  A  specific  technique,  Split-Cycle 


10 


CHAPTER  1.  INTRODUCTION 


Constant  Period  Frequency  Modulation  with  Wing  Bias,  is  introduced  for  tracking  the  wing 
motion  that  provides  a  high  degree  of  control  without  actively  controlling  the  angle  of  attack. 
Particularly  the  frequency  of  the  wing  upstroke  and  downstroke  can  be  varied  to  generate 
a  stroke  averaged  drag  and  a  wing  bias  added  to  provide  pitching  moment  control  all  while 
utilizing  only  a  pair  of  actuators.  Using  the  frequency  modulation  to  control  x  and  z  forcing 
and  the  wing  bias  to  control  moments  and  y  forcing  two  actuators  can  provide  agile  insect 
flight  in  untethered  conditions. 

1.2.2  Flight  Dynamics  of  Flapping  Wing  Systems 

Bidding  [29]  presents  a  flight  dynamic  model  for  flapping  wing  insects  and/or  MAVs.  It 
assumes  a  rigid-body  fuselage  with  wings  that  are  attached  at  a  single  rotating  joint  (that 
move  with  a  prescribed  motion).  The  aerodynamic  model  used  in  his  analysis  is  a  quasi¬ 
steady  model  based  on  the  work  of  Sane  and  Dickinson  [30,  14],  A  spline  interpolation  was 
used  to  approximate  the  observed  data  for  the  wing  kinematics  and  this  data  was  used  to 
generate  the  final  equations  of  motion  (added  to  the  nonlinear  equations  already  derived 
from  the  Newton-Euler  equations).  A  trim  algorithm  was  used  to  calculate  the  periodic 
steady-state  and  Floquet  analysis  used  to  calculate  flight  dynamic  modes  and  stability  to 
small  perturbation.  The  stability  estimate  from  Floquet  analysis  were  close  to  those  from 
the  stroke-averaged  analysis  except  for  slight  differences  in  the  lateral  eigenvalues  and  the 
hovering  case.  This  means  that  the  quasi-steady  approach  may  be  used  for  simplification 
of  the  dynamics  when  necessary.  In  general  the  results  pointed  towards  unstable  or  slightly 
unstable  eigenvalues  necessitating  active  control  but  also  providing  opportunities  for  high 
maneuverability.  Bierling  also  concluded  that  the  close  agreement  between  the  quasi-steady 
analysis  and  the  Floquet  results  might  be  due  to  the  low  wing  mass  of  the  Drosophila  model 
and  that  larger  wing  masses  would  create  larger  differences  in  the  two  analyses.  The  results 
neglect  certain  unsteady  affects  as  part  of  the  quasi-steady  assumption  and  the  results  are 
only  valid  for  the  Drosophila  wing  with  corresponding  flapping  frequency  and  Reynolds 
number. 
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Taylor  et  al.  [31]  examine  the  flight  dynamics  characteristics  of  a  desert  locust.  The 
paper  begins  with  a  brief  discussion  of  flapping  flight  dynamics  in  general  and  concludes 
that  it  is  somewhat  similar  to  forward  rotary  flight  due  to  oscillatory  mass  and  forcing.  The 
authors  discuss  the  distinction  between  open-  and  closed-loop  dynamics  as  pertaining  to 
insects.  Insects  in  general  have  open-loop  dynamics  that  are  not  easily  separated  from  the 
flight  mechanics.  In  general,  blocking  any  kind  of  feedback  usually  results  in  a  stop  in  flight. 
The  authors  go  on  to  discuss  the  characteristics  of  the  desert  locust  including  its  wings  and 
its  flight  sensors.  They  also  provide  analogous  aircraft  sensors  for  many  of  the  locust  sensors, 
for  example,  the  wind  sensitive  hairs  on  the  head  of  the  locust  act  as  angle  of  attack  sensors 
and  side-slip  vanes.  The  Oxford  Zoology  Department’s  low-speed  closed-loop  wind  tunnel 
was  used  to  measure  flight  dynamics  parameters  for  three  different  desert  locusts,  while  two 
others  were  measured  in  an  open-throat  blower  tunnel.  Linear  time-invariant  (LTI)  and 
nonlinear  time-periodic  (NLTP)  models  for  the  flight  dynamics  were  both  analyzed.  The 
models  formulated  by  the  authors  predicted  that  locust  flight  would  be  unstable  without 
visual  feedback  to  the  control  system.  The  NLTP  method  was  also  discussed  as  a  better 
solution  than  the  LTI  model  due  to  the  NTLP  model’s  stronger  forced  oscillatory  character. 
The  authors  also  suggest  that  because  the  MAVs  as  dictated  by  DARPA  will  be  on  the  same 
size  order  as  the  desert  locust,  a  time-periodic  modeling  method  will  be  necessary  to  capture 
their  flight  dynamics  as  well. 

Meirovitch  and  Tuzcu  [32]  formulate  the  equations  of  motion  for  a  flexible  aircraft,  inte¬ 
grating  together  dynamics,  structural  vibrations,  aerodynamics,  and  controls.  Of  particular 
interest  to  this  thesis  was  the  formulation  of  the  dynamic  equations  of  the  aircraft,  treated 
as  a  multi-body  system  using  quasi-coordinates.  The  motions  of  the  body  are  broken  into 
six  degrees  of  freedom  (based  on  the  body  axes  attached  to  the  undeformed  fuselage)  and 
combined  with  deformations  of  the  flexible  components.  Further  detail  on  the  equations  of 
motion  in  terms  of  quasi  coordinates  is  provided  in  Meirovitch  [33]  and  discussed  later  in 
this  thesis. 

Deng  et  al.  [34]  worked  on  mathematical  modeling  of  flapping  flight  MAVs  in  the  inch  size 
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range.  In  this  particular  work  they  worked  on  the  full  system  dynamic  models,  with  particular 
emphasis  on  their  differences  with  respect  to  traditional  fixed  and  rotary  wing  MAVs.  The 
particular  models  included  the  wing-thorax  dynamics,  aerodynamics  (full  flapping  model  at 
low  Reynolds  number),  body  dynamics,  and  a  biomimetic  sensing  system.  The  mathematical 
models  were  derived  from  analytic  solutions,  empirical  models,  and  biological  data.  These 
are  integrated  into  algorithms  for  wing  aerodynamics,  body  dynamics,  actuator,  dynamics, 
external  environment,  and  flight  control,  which  together  make  up  the  Virtual  Insect  Flight 
Simulator  (VIFS).  The  aerodynamics  are  those  of  the  Sane  and  Dickinson  [14]  though  wake 
capture  is  neglected.  The  sensory  system  is  novel  in  that  it  neglects  traditional  flight  sensory 
techniques  in  favor  of  biomimetic  sensors.  Ocelli  are  used  to  estimate  roll  and  pitch  angle, 
while  a  magnetic  compass  estimates  yaw.  Halteres  are  used  for  angular  velocity  sensing  and 
optic  flow  sensors  (much  like  those  described  in  Humbert  et  al.  [35])  are  used  for  navigation 
and  to  avoid  impact.  The  reason  given  for  these  novel  systems  is  that  they  would  be  lighter 
and  draw  less  power  than  conventional  systems,  though  currently  they  are  only  proposed 
and  not  actually  in  use. 

Faruque  and  Humbert  [36]  developed  a  reduced  order  model  for  the  longitudinal  hovering 
dynamics  of  a  fly-like  insect.  They  utilized  a  quasi-steady  aerodynamic  model  first  utilized  by 
Dickinson’s  group  at  Cal  Tech  (averaged  for  each  period  to  a  time  invariant  model)  extended 
by  perturbation  states  from  hover  (to  help  understand  sensing  and  feedback  requirements 
when  moving  away  from  stable  flight)  and  rigid  body  equation  of  motion.  They  then  used  fre¬ 
quency  based  system  identification  to  solve  for  the  transfer  functions  from  the  control  inputs 
to  the  states.  The  heave  dynamics  were  found  to  be  decoupled  from  the  forward  and  pitch 
motions  and  use  of  the  haltere  system  stabilized  the  vehicle  which  otherwise  suffered  from 
an  unstable  oscillatory  mode  (in  agreement  with  previous  CFD  studies).  This  further  rein¬ 
forces  the  idea  of  passive  stabilization  mechanisms  which  help  to  reduce  the  computational 
requirements  of  these  highly  weight  and  power  constrained  vehicles. 

Continuing  work  by  Faruque  and  Humbert  [37]  examined  the  lateral  motion  about  hover. 
The  study  continued  to  examine  passive  stability  mechanisms  and  their  implications  for 
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flight  control  and  stability.  Using  the  previously  described  analysis,  the  lateral  motion  about 
hover  was  shown  to  have  two  stable  poles  and  two  lightly  damped  nearly  unstable  poles. 
One  of  the  poles  (stable)  was  uncoupled  and  corresponded  to  the  yaw  motion  of  the  vehicle 
while  the  others  were  coupled  roll/sideslip  motions.  Again  the  motions  were  found  to  have 
inherent  passive  damping  mechanisms  that  help  to  stabilize  the  vehicle  without  control  input. 
The  designer  can  leverage  these  phenomena  and  reduce  the  burden  placed  on  the  vehicle 
computationally,  reducing  power,  size,  and  weight  issues. 

Sun  and  Xiong  [38]  studied  the  longitudinal  flight  stability  of  hovering  bumblebees  us¬ 
ing  eigenvalue  analysis  to  solve  the  equations  of  motion.  The  aerodynamic  derivatives 
necessary  for  the  analysis  were  computed  using  computational  fluid  dynamics  (CFD).  For 
longitudinally-disturbed  motion  they  identified  three  natural  modes  though  the  instabilities 
arising  from  these  modes  were  all  such  that  the  growing  time  was  much  greater  than  the  wing- 
beat  period  and  as  such  the  bumblebee  would  have  plenty  of  time  to  adjust  wing  motion. 
Sun  and  Xiong  assumed  the  bee  to  be  a  rigid  body  (wings  also  inflexible)  with  six  degrees 
of  freedom  and  the  flapping  wings  were  simplified  by  calculating  the  wingbeat  average  forces 
and  moments  (per  stroke).  The  inertia  effects  of  the  wings  were  also  neglected  (justified  by 
calculating  the  wings’  mass  to  be  about  0.52%  of  the  total  insect  mass).  In  this  particular 
case,  the  equations  of  motion  were  linearized  by  approximating  the  body’s  motion  as  a  series 
of  small  deviations  from  the  reference  condition  (steady,  symmetric  flight).  Sun  and  Xiong 
also  state  that  near  hover  the  body  aerodynamic  forces  are  minuscule  as  compared  to  those 
of  the  wings  and  therefore  can  be  neglected  without  affecting  the  aerodynamic  derivatives 
significantly.  The  eigenvalue  analysis  yielded  one  unstable  oscillatory  mode  and  two  stable 
modes  (subsistence  modes).  Sun  and  Xiong  also  explain  the  physical  meaning  of  these  mo¬ 
tions  and  how  the  bee  could  easily  overcome  the  instability  that  arises  from  the  unstable 
mode. 
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1.2.3  Control  of  Flapping  Wing  Systems 

Hamel  and  Jategaonkar  [39]  present  a  history  of  parameter  identification  in  flight  mechanics. 
It  details  the  development  of  methods  used  in  the  past  as  well  as  currently.  It  also  details 
the  difficulties  in  determining  aircraft  parameters  and  how  those  difficulties  have  changed 
over  the  course  of  time.  Some  of  the  early  methods  presented  were  the  first  attempts  as 
determining  aerodynamic  forces  and  included  dynamic  flight  tests,  the  Longitudinal  Oscilla¬ 
tion  Method,  the  Pulse  Method  of  Dynamic  Response  Testing,  the  Time-Vector  Method  and 
Analog  Matching.  Most  of  these  methods  have  become  outdated,  even  at  the  time  of  the 
publication  of  this  paper,  and  Hamel  and  Jategaonkar  go  on  to  present  more  modern  meth¬ 
ods  of  parameter  estimation.  Most  of  these  methods  were  brought  about  with  the  advent 
and  increase  in  power  of  computers,  especially  the  personal  computer.  They  also  discuss 
four  important  aspects  of  parameter  identification  in  modern  methods,  namely:  design  of 
the  control  input  shape,  selection  of  instrumentation  and  filters  for  high  accuracy,  definition 
of  the  structure/mathematical  model,  and  proper  selection  of  the  time  or  frequency  domain 
identification  method.  Each  of  these  areas  are  discussed  in  depth  and  specific  examples  given. 
There  is  also  a  discussion  on  system  identification  applied  to  related  topics.  The  paper  goes 
on  to  conclude  that  the  current  techniques  presented  in  the  paper  provided  sufficiently  ac¬ 
curate  results  and  had  been  tested  thoroughly  enough  to  be  considered  powerful  tools  not 
only  for  research  but  also  for  industry.  It  also  hypothesizes  that  the  area  of  parameter  iden¬ 
tification  will  continue  on  limited  only  by  the  imagination  and  innovation  of  analysts,  and 
motivated  by  the  need  to  better  understand  aerodynamic  phenomena. 

Humbert  et  al.  [35]  deals  with  the  nature  of  insect  sensory  feedback  control  and  how  it 
could  be  implemented  in  MAVs.  The  primary  mechanism  most  desirable  in  insect  sensory 
motor  interaction  is  that  of  sensorimotor  convergence.  The  output  of  sensory  neurons  is 
conducted  directly  to  motor  control  centers  rather  than  central  processing  areas  for  faster 
feedback  than  possible  by  first  processing  and  then  sending  to  motor  centers.  The  visual 
sensory  output  focused  on  by  this  paper  is  optic  flow,  which  is  essentially  the  derivative  of  the 
objects  in  the  optic  field.  The  insect’s  eyes  process  the  velocity  of  the  objects  in  its  optic  field 
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and  also  their  proximity.  This  minimized  data  is  passed  to  the  motor  control  centers  for  quick 
reactions  and  maneuvering.  This  is  what  gives  insects  their  characteristic  maneuverability 
and  quick  reaction  times.  Ideally,  this  data  would  also  include  the  spatial  structure  and 
distribution  of  objects  in  the  flow  field  to  aid  in  navigation.  Humbert  et  al.  go  on  to 
describe  the  process  of  Wide-Field  Integration  (WFI)  processing  of  ideal  planar  optic  flow, 
including  characterization  and  interpretation  of  WFI  outputs.  They  also  describe  how  output 
feedback  of  WFI  sensory  information  can  be  used  to  create  pitch  altitude  stability  and  terrain 
following  capability  (demonstrated  using  planar  rotorcraft  dynamics).  It  is  concluded  that 
terrain  following  applications  to  date  have  only  utilized  a  small  amount  of  the  information 
available  from  optic  flow.  It  would  be  possible  to  significantly  improve  terrain  following, 
closed  loop  stability,  and  performance  through  proper  implication  of  optic  flow  sensing.  They 
demonstrate  a  methodology  of  stability  through  optical  flow  data  with  a  minimum  of  sensory 
information  from  other  methods.  For  instance  only  a  measurement  of  forward  speed  was 
necessary  to  achieve  zero  steady  state  error  in  tracking  altitude  reference,  and  only  pitch 
orientation  necessary  for  hover  equilibrium. 

Hedrick  et  al.  [40]  discuss  passive  yaw  stabilization  in  a  wide  range  of  flapping  fliers. 
They  propose  that  rather  than  actively  damping  out  oscillations  and  perturbations,  flying 
organisms  use  passive  stability  which  damp  out  angular  velocity  through  a  methodology 
the  authors  term  flapping  counter  torque  (FCT).  Their  model  predicts  similar  damping  on 
a  per-wingbeat  time  scale  for  isometrically  scaled  animals.  They  show  that  animals  may 
specialize  in  both  maneuverability  and  stability  despite  the  normal  assumption  that  these 
two  quantities  run  counter  to  each  other.  They  also  propose  that  other  mechanisms  besides 
yaw  could  be  explained  by  this  phenomenon.  The  particular  turn  studied  in  this  paper 
was  a  low-speed  yaw  turn  of  greater  than  or  equal  to  60  degrees.  This  turn  was  described 
as  having  an  angular  acceleration  phase  (with  active  control)  and  an  angular  deceleration 
phase(with  either  passive  or  active  control).  Traditionally  larger  animals  such  as  birds  are 
modeled  as  requiring  active  deceleration  to  overcome  their  larger  inertia  due  to  mass  while 
smaller  animals  such  as  insects  are  modeled  with  significant  fluid  drag,  allowing  for  passive 
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deceleration.  Considering  an  animal  in  hover  or  low  speed  flight,  damping  from  wing  motion 
arises  from  drag  created  by  the  flapping  motion.  When  the  animal  is  in  full  body  rotation 
(yaw),  net  outside  wing  velocity  is  increased  on  the  downstroke  and  net  inside  wing  velocity 
increased  on  the  upstroke.  This  gives  rise  to  a  torque  that  slows  the  animal’s  rotation 
(on  both  strokes  the  net  torque  due  to  force  imbalance  acts  to  slow  the  rotation).  The 
authors  refer  to  this  torque  as  FCT  because  it  is  flapping  induced  and  acts  counter  to  the 
already  present  rotation.  Hedrick  et  al.  then  go  on  to  develop  an  equation  to  model  FCT 
and  another  equation  for  active  deceleration  through  torque  generation.  They  use  the  two 
equations  to  predict  rotational  deceleration  dynamics  and  compare  these  predictions  to  actual 
measurements  of  yaw  turning  in  seven  species  of  flying  animals  over  six  orders  of  magnitude  of 
mass.  The  FCT  model  was  able  to  accurately  predict  yaw  deceleration  over  all  seven  species 
(four  invertebrate  and  three  vertebrate)  while  the  active  model  was  not  able  to  accurately 
predict  the  deceleration  rate  measured  in  these  species.  While  active  control  is  possible,  the 
authors  contend  that  the  primary  mechanism  in  yaw  turn  deceleration  is  passive.  Due  to  these 
results,  maneuverability  and  stability,  often  cast  in  opposition  to  each  other,  are  more  readily 
optimized  together.  Essentially,  enhancing  some  factors  that  enhance  maneuverability  would 
also  increase  FCT,  itself  an  active  form  of  passive  stability.  One  such  factor  discussed  in  this 
paper  is  wingbeat  frequency.  Increasing  frequency  would  allow  for  greater  maneuverability 
and  active  stabilization  as  well  as  increasing  FCT  through  greater  force  production.  This 
does,  however,  come  with  the  cost  of  increases  in  power  requirements  to  achieve  these  higher 
flapping  frequencies.  The  authors  also  note  that  while  FCT  does  provide  open  loop  stability, 
reducing  neuromuscular  and  sensory  requirements,  it  does  not  eliminate  these  as  FCT  results 
in  asymmetric  forcing  in  the  wings  which  suggests  neurological  input  to  achieve  symmetry 
of  the  flapping  motion.  The  authors  also  discuss  the  fact  that  fast  forward  flight  is  also 
likely  influenced  by  FCT  but  pitching  and  longitudinal  dynamics  are  known  to  be  inherently 
unstable  in  open  loop  conditions  and  require  active  control. 

Doman  et  al.  [41]  developed  a  method  of  controlling  a  flapping  wing  micro  air  vehicle 
through  varying  the  velocity  profiles  of  each  wing.  The  rotation  of  the  wing  is  a  passive  mech- 
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anism  while  the  stroke  motion  is  actuated  using  what  is  termed  a  Split-Cycle  Constant-Period 
Frequency  Modulation.  This  was  accomplished  through  the  use  of  a  split  cycle  parameter 
that  shifted  the  peak  of  the  flapping  cycle  while  maintaining  the  same  period.  Through  the¬ 
oretical  analysis  and  simulation  this  method  of  control  is  shown  the  be  capable  of  of  control 
over  the  vertical  and  horizontal  body  forces  as  well  s  the  rolling  and  yawing  moment  when 
just  controlling  each  wing’s  flapping  frequency  and  the  split  cycle  parameter.  By  introducing 
a  bob-weight  pitching  moment  can  also  be  controlled.  The  controller  is  implemented  using 
derived  sensitivities  of  cycle-averaged  forces  and  moments  (from  changes  in  the  wing  beat 
kinematic  parameters).  These  sensitivities  are  used  to  formulate  a  cycle-averaged  control  law 
that  successfully  stabilizes  two  different  models,  one  using  blade-element  theory  and  another 
based  on  an  empirical  unsteady  aerodynamic  model  derived  from  experiments. 

Liu  et  al.  [42]  developed  a  computational  framework  to  provide  future  guidelines  for  MAV 
design.  The  framework  integrated  aerodynamics,  flight  dynamics,  structural  dynamics,  and 
flight  stability  and  maneuverability  applied  to  hawkmoth  and  fruit  fly  models.  The  aero¬ 
dynamic  model  was  a  fully  unsteady  Navier-Stokes  grid  solver  coupled  to  a  finite  element 
solver  to  account  for  the  large  deformations  of  flexible  wings  (results  provided  for  both  rigid 
and  flexible  wings).  The  flight  dynamic  model  is  a  six-degree-of-freedom  set  of  dynamic 
equations  (Newton-Euler  scheme)  with  translation  described  in  the  inertial  frame  and  rota¬ 
tion  described  in  the  body  frame.  In  the  end,  all  of  these  are  coupled  through  a  non-linear 
coupling  model  and  simulated  in  order  to  aid  in  analysis  of  bio-inspired  flight. 

Continuing  work  from  [34],  Deng  et  al.  [43]  present  research  on  flight  control  algorithms 
inspired  by  the  top-down  architecture  of  real  insects  in  order  to  try  and  preserve  high  per¬ 
formance  aspects  despite  limited  computing  power.  The  hierarchical  nature  of  the  model 
allows  for  modular  construction  of  the  controller,  with  each  module  responsible  for  its  own 
particular  task.  The  independent  nature  of  each  module  allows  for  later  improvement  with¬ 
out  redesign  of  the  entire  system.  The  main  levels  considered  are  the  navigation  planner,  the 
flight  mode  stabilizer  and  the  wing  trajectory  controller.  The  navigation  planner  uses  data 
from  sensors  and  commands  to  choose  the  sequence  of  flight  modes  necessary  to  achieve  the 
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mission  and  navigate  the  terrain  around  the  vehicle.  These  flight  modes  are  passed  to  the 
flight  mode  stabilizer  which  uses  inputs  from  the  biomimetic  sensors  described  in  Deng  [34]  to 
stabilize  the  various  flight  modes.  It  then  chooses  the  necessary  forces  and  torques  to  achieve 
this  stable  configuration  and  sends  them  to  the  wing  trajectory  controller  which  chooses  the 
appropriate  input  voltages  for  the  actuators.  The  particular  advances  in  this  work,  at  that 
time,  came  in  the  form  of  the  approximation  of  the  time  varying  aerodynamics  as  a  time 
invariant  system  using  averaging  theory  (through  assuming  the  wing-beat  frequency  is  much 
higher  than  that  of  the  body  motion,  meaning  that  only  the  mean  forces  and  moments  over 
one  aerodynamic  period  would  significantly  affect  the  body  dynamics)  and  in  the  use  of  a 
biomimetic  parametrization  (such  as  timing  of  rotation,  stroke  angle  amplitude,  stroke  angle 
offset,  etc)  of  the  individual  wing  paths.  These  simplifications  allow  the  controller  to  take 
take  data  from  on  board  sensors  and  voltage  inputs  to  the  wing  actuators  and  control  the 
system,  without  the  complication  of  the  controller  needing  to  know  the  exact  aerodynamic 
model  or  accurate  vehicle  morphological  parameters. 

Oppenheimer  et  al.  [44]  use  their  work  from  [28]  and  develop  a  control  strategy  to 
allow  the  vehicle  to  hover  stably  using  two  wing  actuators  controlling  the  angular  position 
of  the  wing  in  the  stroke  plane  and  a  bob  weight  to  move  the  vehicle  center  of  gravity. 
The  Split-Cycle  Constant-Period  Frequency  Modulation  technique  is  used  to  generate  non 
zero  stroke  averaged  rolling  and  yawing  moments  while  the  bob  weight  generates  pitching 
moments  and  controls  translation.  The  bob  weight  had  to  be  included  as  using  the  Split- 
Cycle  Constant-Period  Frequency  Modulation  technique  to  generate  pitching  moments  also 
generated  vertical  forcing  and  rolling  moments  whereas  the  bob  weight  could  generate  the 
pitching  moment  independent  of  other  forcing.  Expressions  for  the  control  derivatives  due 
to  variations  in  the  wing  beat  frequencies,  split  cycle  parameters,  and  movement  of  the 
bob  weight  are  given.  Using  the  stroke-averaged  aerodynamic  model,  a  flight  control  law  is 
generated  and  then  applied  to  the  time  varying  dynamic  model  to  demonstrate  applicability 
of  the  stroke  averaged  model. 
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Equations  of  Motion 


In  this  chapter  the  equations  of  motion  are  derived  using  quasi  coordinate  form  of  the  typical 
Lagrange’s  equations.  The  purpose  and  uses  of  quasi  coordinates  are  described  along  with 
how  they  are  applied  to  this  particular  system.  The  changes  to  the  Lagrange’s  equation  due 
to  the  use  of  quasi  coordinates  rather  than  the  usual  general  coordinates  are  also  presented 
with  more  detail  available  in  Meirovitch  [33].  This  is  followed  by  a  formulation  of  the  kinetic 
and  potential  energy  for  use  in  the  modified  Lagrange’s  equations  (along  with  the  virtual 
work  due  to  the  non-conservative  forces).  The  Lagrange’s  equations  are  then  fully  formulated 
and  transformed  into  a  state  space  form  for  ease  of  use.  The  wing  kinematics  used  in  this 
system  are  also  described,  followed  by  a  check  of  the  linear  and  angular  momentum  due  to 
only  the  inertial  loading.  The  actual  physical  quantities  used  to  describe  this  system  are 
provided  in  Section  5.1. 

2.1  Quasi  Coordinates 

Quasi  coordinates  are  a  set  of  independent  coordinates  that  describe  a  system  uniquely.  In 
the  case  of  a  rigid  system  these  coordinates  are  often  those  that  describe  the  translation  of 
the  body  origin  and  the  rotation  of  the  reference  axes  [33].  In  the  case  of  this  thesis  these 
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coordinates  are  the  velocity  vector  and  position  vector  of  the  body  frame,  b ,  relative  to  an 
inertial  frame,  i,  and  the  angular  velocities  and  angles  of  the  body  frame  relative  to  the 
inertial  frame  (as  described  below  in  Equation  2.1  through  Equation  2.4).  These  coordinates 
are  described  as  quasi  coordinates  because  the  velocities  can  not  be  directly  integrated  to 
obtain  positions/angles. 

Figure  2.1  illustrates  the  system  investigated  in  this  thesis.  The  system  is  described  by 
several  frames  including  the  inertial  frame  (center  at  origin,  0),  i,  the  body  frame  (center  at 
body  center  of  mass,  B),  b ,  the  stroke  frame,  s,  and  the  wing  frame  (center  at  wing  center 
of  mass,  W),  w.  Using  these  frames  the  quasi  coordinates  can  be  described  as  follows: 

{lrB}  =  {lxBnyB,lzB}T  (2.1) 

{<pb}  =  {T6,  06,  (2.2) 
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{6  vB}  =  {buB  ,bvB  ,bwB}T  (2.3) 

{bub}  =  {bpb,bqb,brb}T  (2.4) 


Here  the  left  subscript  denotes  the  reference  frame  the  variable  is  written  in  while  the  upper 
right  superscript  denotes  the  point (s)  or  frame  being  referenced. 


In  order  to  relate  the  velocities  and  angular  velocities  to  the  position  and  angles  one  must 
first  take  into  account  the  rotation  between  the  body  and  inertial  frame  through  the  Euler 
angles  described  in  Equation  2.2.  This  is  a  simple  3-2-1  rotation  giving  rise  to  the  rotation 
matrix: 


[Tk] 


cos  T6  cos  Qb 

cos  T6  sin  @b  sin  —  sin  T  cos  <E>b 
cos  T6  sin  06  cos  +  sin  T6  sin  ‘h6 


sin  T6  cos  0b 

sin  T6  sin  Qb  sin  <f>b  +  cos  T6  cos 
sin  T6  sin  Qb  cos  <3>b  —  cos  T6  sin  <f>6 


—  sin  Qb 

cos  06  sin  <f>6 

cos  0 b  cos 
(2.5) 


and  the  position  and  angles  are  then  related  to  the  velocities  and  angular  velocities  by: 


{A5}  =  [Tib]{fevB}  (2.6) 

=  [E  «]-1W}  (2.7) 


where  Ebi  is  the  matrix  that  relates  the  body  angular  rates  and  the  time  derivatives  of  the 
Euler  angles.  The  derivation  of  this  matrix  can  be  found  in  the  appendix  to  Bierling  [1], 


[Eu] 


—  sin  0 
sin  $  cos  0 
cos  $  cos  0 


0  1 

cos  $  0 

—  sin  0 


(2.8) 


For  convenience  a  new  plane  is  defined,  the  stroke  plane,  in  which  the  majority  of  the 
flapping  motion  takes  place.  This  plane  is  explained  by  Ellington  [27]  as  the  plane  defined 
by  the  wing  bases,  that  is,  the  root  of  the  wing  and  the  wing  tip  at  its  forward  and  rearmost 
points.  This  definition  also  allows  for  a  differentiation  of  the  stroke  plane  for  each  wing 
which  could  be  necessary  for  control  and  maneuvering  as  shown  by  Fry  [6].  The  stroke  plane 


22 


CHAPTER  2.  EQUATIONS  OF  MOTION 


Figure  2.2:  Rotation  from  the  stroke  frame  to  the  wing  frame[l] 


is  defined  by  a  2  rotation  about  the  body  y  axis  through  the  angle  Sr  or  Si  (depending  on 
whether  it  is  the  right  or  left  wing).  Thus  the  rotation  matrix  from  the  body  to  the  stroke 
frame  is: 


[T,J 


cos  S  0  —  sin  5 
0  1  0 

sin  5  0  cos  S 


(2.9) 


The  wing  frame  is  then  defined  relative  to  the  stroke  frame  by  three  Euler  angles,  namely 
the  stroke  angle  k,  the  rotation  angle  r,  and  the  deviation  angle  /i.  The  rotation  itself  is  a 
3-1-2  rotation  through  n  then  ji  and  finally  r  as  shown  in  Figure  2.2.  The  rotation  matrix 
from  the  stroke  frame  to  the  wing  frame  is  therefore: 


—  sin  k  sin  /i  sin  r  +  cos  k  cos  r  cos  n  sin  /i  sin  r  +  sin  k  cos  r  —  cos  y  sin  r 
—  sin  k  cos  y  cos  n  cos  //  sin  /i 

sin  k  sin  /i  cos  r  +  cos  n  sin  r  —  cos  n  sin  /i  cos  r  +  sin  n  sin  r  cos  y  cos  r 

(2.10) 


The  angular  velocity  of  the  wing  frame  and  relative  to  the  body  frame  is  given  by  the 
sum  of  the  angular  velocity  of  the  stroke  frame  relative  to  the  body  frame  and  the  angular 
velocity  of  the  wing  frame  relative  to  the  stroke  frame  (details  of  this  derivation  are  given  by 
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Bierling  [1]): 

{ubw}  =  {cjbs}  +  {w™}  (2.11) 

The  angular  velocity  of  the  stroke  frame  relative  to  the  body  frame  is  simply: 

{VS}  =  {0,<5,0}T  (2.12) 

while  angular  velocity  of  the  wing  frame  relative  to  the  stroke  frame  is  given  in  both  frames 
by: 


>  (2.13) 


f  y 

— f  sin  k  cos  fi  +  fi  cos  k 

'  ' 

—k  sin  r  cos  fi  +  fi  cos  r 

{sV™}  =  < 

f  cos  n  cos  fi  +  fi  sin  n 

>  {w™}  =  < 

k  sin  fi  +  f 

f  sin  a  +  k 

V  '  > 

k  cos  r  cos  u  +  u  sin  r 

Taking  the  time  derivative  yields  the  following  vector  which  will  be  needed  later  in  the 
derivation  of  the  equations  of  motion  and  the  aerodynamic  model. 


on  = 


— k  sin  r  cos  fi  —  kf  cos  r  cos  fi  +  kfi  sin  r  sin  fi  +  jl  cos  r  —  /if  sin  r 

k  sin  fi  +  kji  cos  fi  +  f 

k  cos  r  cos  fi  —  kf  sin  r  cos  fi  —  kfi  cos  r  sin  fi  +  ji  sin  r  +  /if  cos  r 


>  (2.14) 


2.1.1  Quaternions 


Quaternions  represent  an  alternative  to  the  traditional  Euler  angle  approach  for  equations  of 
motion.  Utilizing  four  parameters  instead  of  three,  quaternions  avoid  the  gimbal  lock  problem 
that  Euler  angles  exhibit  as  the  pitch  angle  approaches  ninety  degrees.  However,  quaternions 
also  exhibit  the  disadvantage  of  providing  less  physical  insight  due  to  their  abstract  nature. 
In  this  study  quaternions  are  defined  as  follows: 


Q  =  <7o  +  qii  +  Q2j  +  qzk  (2.15) 

The  rotation  matrix  from  the  body  frame  to  the  inertial  frame  (and  visa  versa)  as  well  as 
the  kinematic  equations  can  be  written  in  terms  of  quaternions  as  described  in  [45]. 
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Tbi  — 


9o  +  9?  ~  Q2  ~  ?3 

2  (qiqz  ~  9o93) 

2  (qm  +  q 092) 


2  (91 92  +  9o?3) 
9o  -  9?  +  q\  -  q\ 

2  (9293  -  9o9i) 


2  (9193  -  90^2) 

2  (92^3  +  9o9i) 

ql  -  q\  -  q\  +  ql 


(2.16) 


-q  1 

-92 

-93 

<?0 

-93 

92 

<?3 

9o 

— 9i 

-92 

9i 

9o 

(2.17) 


2.2  Lagrange’s  Equations  for  Quasi  Coordinates 


The  derivation  of  the  equations  of  motion  is  based  on  the  standard  Lagrange’s  equations. 
These  are  usually  derived  from  Hamilton’s  principle  and  are  valid  for  a  set  of  generalized  co¬ 
ordinates.  As  this  is  a  well  established  derivation,  it  will  not  be  repeated  here  but  more  detail 
is  provided  in  Meirovitch  [33].  Typically  the  standard  Lagrange’s  equations  are  expressed  as: 


d  f  dL\  dL 

M\Wk)~Wk=Qk 


k  =  1,2 


(2.18) 


where  L  is  the  Lagrangian  of  the  system  (difference  of  the  kinetic  and  potential  energy), 
qk  are  the  generalized  coordinates,  and  Qk  are  the  generalized  (typically  non-conservative) 
forces. 


The  system  investigated  in  this  thesis  is  described,  however,  by  a  set  of  quasi  coordinates, 
necessitating  that  a  different  form  of  Lagrange’s  equations  be  used.  Lagrange’s  equations 
as  applied  to  quasi  coordinates  are  derived  by  Meirovitch  [33].  The  general  form  of  these 
equations  is  as  follows: 


d 

dt 


=  {N} 


(2,19) 


where  q  is  the  generalized  coordinate,  u  is  a  quasi  coordinate  representing  of  the  velocities 
(cannot  be  directly  integrated  to  get  the  positions  q).  and  L  is  the  Lagrangian  written  in  a 
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new  form  as  a  function  of  the  coordinates  qb  and  Uk-  Because  the  u's  can  be  written  as  linear 
combinations  of  the  q\ .  terms  they  can  be  expressed  in  the  form  ujs  =  aisqi+ct2sq2  +  ■■■  +  Oinsqn 
or: 

{cu}  =  [a]T{q},  {q}  =  [/3]{cu}  (for non  —  singular  [a])  (2.20) 


Expressing  equation  2.19  in  terms  of  the  variables  relevant  to  this  system  yields: 
d  f  dL  \  ,  _  ,,,  dL  r_  „  dL 


d 


(. 


dL 


dt  \d{bvB} 


+  {&b} 


d{bvB} 


■  bi\ 


d{iYB} 


=  {b?B} 


dL 


+  {&b}- 


dL 


fE 


dL 


dt\d{bub}J  '  l°~Jd{bvB}  '  l°~Jd{b u:b}  l^bild{i(pB} 
where  L  is  the  Lagrangian  of  the  system  and  bFB  and  foMBare  the  aerodynamic  and  grav¬ 
itational  forces  and  moments  respectively  (expressed  in  the  body  frame  b ).  Because  all  of 
the  bodies  in  this  system  are  rigid  and  the  gravitational  effects  are  handled  as  external  loads 
there  is  no  potential  energy  involved  in  the  Lagrangian  calculation.  Also,  because  the  kinetic 
energy  is  not  dependent  on  the  positions  {irB}  or  the  angles  {i<pB}  the  corresponding  terms 
in  equations  2.21  and  2.22  are  both  zero.  More  details  on  the  exact  form  of  the  kinetic  energy 
equation  follow. 


=  {Mb} 


(2.21) 


(2.22) 


2.2.1  Kinetic  Energy  Equation 


The  expression  for  the  kinetic  energy  is  typically  given  as  one  half  a  mass  component  mul¬ 
tiplied  by  velocity  squared.  Expanding  this  into  a  little  more  robust  form  provides  the 
equation: 


KE 


1 

2 


\T\ dm 


(2.23) 


and  for  this  particular  system  (the  equations  for  the  second  wing  are  derived  in  the  same 
exact  manner  as  for  the  first): 


KE  =  - 


VB  V1 


dmb  +  - 


V'v  V"  dm „ 


(2.24) 


After  integration  this  can  be  condensed  to  a  convenient  matrix  equation: 


KE  =  1{Vt}[M]{V} 


(2.25) 
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In  order  to  preserve  this  form  the  vector  V  must  be  enlarged  to  include  the  rotational  velocity 
between  the  body  and  wing  frames.  This  vector  is  given  by: 

{V}  =  {bvBT  ,bubT  ,bubwT}T  (2.26) 


This  allows  us  to  create  a  matrix  M  as  defined  below: 


[M] 


Mn 

Mi2 

Mi3 

M2i 

m22 

m23 

M3i 

m32 

m33 

[Mn]  = 

%[  E 

where  [E]  is  the  3x3  identity  matrix. 


(2.27) 


(2.28) 


[M12]  =  mw  (-brBH  -  TbwwrHWTwb)  (2.29) 

[M13]  =  mw  (- TbwwrHWTwb )  (2.30) 

[M21]  =  [M12]t  (2.31) 

[M22]  =6  if  +  TbwwI^Twb  +  mw  (- brBH  -  Tbwwr™  Twb)T  ( -brUH  -  Tbwwr1IU  Twb ) 

(2.32) 

[M23]  =  TbwwI™Twb  +  mw  (-brBH  -  TbwwrHWTwb)T  (- TbwwrHWTwb )  (2.33) 

[M31]  =  [M13]t  (2.34) 

[M32]  =  [M23]t  (2.35) 

[M33]  =  TbwwI^Twb  +  mw  (~TbwwrHWTwb)  (—Tbwwrhn  Twft)  (2.36) 

The  derivatives  with  respect  to  the  velocities  v6  and  u:b  are  as  follows: 


^  —  -  {bVB  Mn  +6  M2i  +b  cjbu  M3i  j  +  -  (MnbVB  +  Mi2(,aji  +  Mi^ii/"’)  (2.37) 

or  -j  i 

—  2  ( V Mi2  +b  u)b  M22  +b  u)bw  M32)  +  -  (M2ifeVB  +  M225o;6  +  NL23bujbw)  (2.38) 
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and  their  time  derivatives: 


Jt  ( =  2  B'Ml 1  +b  ^  M21  +fc  ^  M21  +fc  ^  Msi  +fe  “>}W'  M 

+  1  (M116vB  +  M12bub  +  Ml2bub  +  M13bu;bw  +  M13bubw 
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(2.39) 


d  f  dL  \ 

Jt  \JJJ) 


=  ^  (6VfiTM12  +fe  vsTM12  +6  ubTM22  +b  u;bTM22  +b  u>bwT M32  +b  d^TM32) 

+  -  ^M2ihVii  +  M2if,v'B  +  Jl22bujb  +  'Wl22bujb  +  M23fcu;6ti’  +  ~M.23b&bw j 

(2.40) 


With  these  equations  it  is  possible  to  begin  assembling  the  full  form  of  the  Lagrange’s  equa¬ 
tions  stated  previously  in  Equations  2.21  and  2.22.  It  is  important  to  note  that  the  equations 
of  motion  also  require  the  external  forces  and  moments.  For  the  purpose  of  this  thesis  these 
would  include  the  gravitational  effects  and  those  effects  due  to  the  aerodynamics  of  the  wings 
(the  body  is  not  considered  to  have  an  aerodynamic  affect). 


2.2.2  State  Space  Form 

For  ease  of  use  the  equations  of  motion  are  transformed  into  a  state  space  form  in  which  the 
body  acceleration  and  body  angular  acceleration  terms  are  collected  on  the  left  hand  side 
of  the  equation  while  all  others  are  kept  on  the  right  hand  side.  This  yields  the  following 
equations: 

O  TT 

MnvB  +  M12ih6  =  Fa  +  Fg-  ub—  -  (M12u>&  +  M13ub™  +  M13u;to)  (2.41) 

d  L  dL  f  \ 

M2ivb  +  M22ci;6  =  +  Mg  —  wB  — — —  ujb  — — ^  —  ^M2iWB  +  ~Wl22ujb  +  M23u^1"  +  M23u^u  j 

(2.42) 

In  simplified  terms  these  equations  are  of  the  form  of  a  mass  matrix  dependent  on  the 
kinematic  variables  (which  are  periodic  with  time)  multiplied  by  the  state  derivatives  equal 
to  the  forcing  function  dependent  on  the  states  and  the  kinematic  variables: 
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[M(C(t))]x  =  f(x,C(t)) 

Where  x  are  the  states  and  £(t)  are  the  periodic  prescribed  kinematics. 


(2.43) 


2.2.3  Center  of  Gravity  Calculation 

The  center  of  gravity  of  the  entire  system  is  not  centered  about  the  body  but  rather  moves 
with  respect  to  the  body  due  to  the  motion  of  the  wing.  The  position  of  the  center  of  gravity 
with  respect  to  the  inertial  frame  is  calculated  as  (following  equations  only  take  into  account 
one  wing,  the  second  wing  is  added  in  the  same  fashion): 


B  i  W 

CG  iV  mB  +i  r  mw 

jT*  =  - 

rri  R  +  mw 

Expanding  trM  into  its  various  components  yields: 


(2.44) 


,CG 


iVBmB  +  mw  {iTB  +  Tib  brBH  +  TibTbw  wrmi ) 
rriT 


(2,45) 


The  kinematic  equation  for  the  motion  of  the  center  of  gravity  is  simply  the  derivative  of 
Equation  2.45. 


,r 


CG 


TibbvBmB  +  mw  ( TibbvB  +  Tib  bub  brBH  +  Tib  bL jbTbw  wrHW  +  Tib  bu >bwTbw  wrHV[ ) 

rriT 

(2.46) 


2.2.4  Wing  Kinematics 

The  wing  kinematics  are  divided  into  three  main  motions,  the  dominant  motion  being  the 
stroke  (often  referred  to  in  terms  of  the  downstroke  and  upstroke),  the  rotation  of  the  wing 
about  its  axis  of  rotation  (to  ensure  that  the  wing  is  at  a  positive  angle  of  attack  for  all  or  the 
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majority  of  the  stroke),  and  the  deviation  from  the  stroke  plane.  As  insects  do  not  employ 
a  tail  for  steering  or  stabilization,  the  wing  kinematics  are  of  particular  importance  since 
they  are  responsible  for  the  aerodynamic  forcing,  and  passive  as  well  as  active  stabilization. 
The  focus  of  this  thesis  was  not  to  explore  the  exact  kinematics  of  a  particular  insect  and 
therefore  a  more  simple  approach  was  used  to  describe  the  wing  motion  and  help  with  physical 
understanding  of  the  phenomena  occurring  during  flight.  In  this  case  a  sinusoidal  description 
of  the  wing  motion  determined  through  use  of  an  amplitude,  offset,  and  phase  lag  captures 
the  necessary  motion  well  enough. 

2.2.4. 1  Stroke  Angle 

The  stroke  angle  is  the  dominant  motion  of  the  wing  which  for  many  insects  is  a  rowing 
motion  in  hover  rather  than  the  typical  up  and  downstroke  of  birds  (though  some  insects 
such  as  some  butterflies  and  moths  utilize  a  motion  dominated  by  up  and  down  strokes).  In 
this  case  the  motion  is  described  by  the  following  equation  (where  v  is  the  flapping  frequency): 


K  =  hio  +  Umax  COs(27 TUt)  (2.47) 

These  variables  (k0,  Kmax)  are  separate  for  each  wing  where  setting  Kqx  =  — ^0,z  and 
Kmax,r  —  —Kmax,i  produces  a  symmetric  stroke  pattern,  an  example  of  which  is  shown  below 
in  Figure  2.3. 

2. 2. 4. 2  Deviation  Angle 

Choosing  a  frequency  twice  that  of  the  stroke  and  rotation  provides  for  the  figure  eight 
pattern  characteristic  of  many  insect  flapping  motions.  The  figure  eight  pattern  is  illustrated 
in  Figure  2.4  where  it  is  helpful  to  note  that  the  stroke  plane  is  the  zero  line  (in  z)  and  all 
motion  above  or  below  this  line  is  due  to  deviation.  In  order  to  produce  this  motion  it  is 
necessary  to  introduce  a  phase  lag  in  deviation,  relative  to  the  stroke  angle. 
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Figure  2.3:  Example  stroke  angle 


x  10"3 


Figure  2.4:  Right  wing  tip  location  in  the  x-z  plane,  showing  characteristic  figure  eight  pattern 
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Figure  2.5:  Example  deviation  angle 


fi  =  /_i0-  /_imax  cos(47tzA  +  (2.48) 

Again  this  motion  is  separate  for  each  wing  and  setting  /r0>r  =  —  /r0,z,  4>^,r  =  and 
Hmax,r  =  ~^max,i  produces  a,  symmetric  stroke  pattern,  an  example  of  which  is  shown  in 
Figure  2.5. 

2. 2. 4. 3  Rotation  Angle 

The  rotation  angle  produces  much  of  what  would  typically  be  considered  the  geometric  angle 
of  attack  and  is  also  lagged  from  the  stroke  angle,  this  time  by  <j>T.  This  lag  allows  for 
advanced  or  delayed  rotation,  where  a  negative  phase  lag  produces  an  advanced  rotation 
(rotation  begins  before  the  stroke  reaches  the  maximum  amplitude  at  the  end  of  each  half 
stroke). 


r  =  r0  -  rmax  sin(27 Tvt  +  0T)  (2.49) 

Separating  the  motion  for  each  wing  involves  having  different  TotTmaXi  and  <f>T  where 
To,r  =  ro ,i  <Ar,r  =  and  Tmax,r  =  Tmaxj  produces  a  symmetric  stroke  pattern,  an  example  of 
which  is  shown  in  Figure  2.6. 
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Figure  2.6:  Example  rotation  angle 

2.2.5  Conservation  of  Momentum 

Applying  the  wing  kinematics  discussed  in  section  2.2.4  in  a  vacuum  without  gravity  creates 
an  environment  where  no  external  loads  are  considered.  This  allows  the  consideration  of  the 
system  with  only  the  inertial  forcing  from  the  wing  motion  and  can  be  used  for  validation 
the  dynamics  in  vacuum. 

The  initial  condition  determines  the  initial  momentum  which  is  conserved  because  there 
is  no  external  forcing.  Setting  the  initial  state  (vB  and  ub)  equal  to  zero  will  give  a  non-zero 
initial  linear  and  angular  momentum  so  it  is  more  convenient  to  first  solve  for  an  initial 
state  whereby  the  linear  and  angular  momentum  are  equal  to  zero.  The  linear  and  angular 
momentum  equations  below  are  the  result  of  the  principles  of  linear  and  angular  momentum, 
for  more  details  see  Bierling  [1], 

LM  =  vB  (mb  +  mw)  +  mw  (ub  x  rBW  +  ubw  x  rHW )  (2.50) 

AM  =  Ifw6  +  l{£  (ub  +  ubw)  +  rBW  x  (vB  +  ub  x  rBW  +  x  rHW )  mw  (2.51) 

Setting  these  equations  equal  to  zero  and  solving  for  the  resulting  initial  state  (vB  and  ujb) 
gives  a  point  from  which  time  marching  can  begin.  Simulating  the  system  and  calculating  the 
resulting  linear  and  angular  momentum  (in  the  inertial  frame)  provides  Figures  2.7  and  2.8. 
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Figure  2.7:  Linear  momentum  in  vacuum 
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Figure  2.8:  Angular  momentum  in  vacuum 


This  shows  that,  as  expected,  the  linear  and  angular  momentum  are  equal  to  zero  withing 
numerical  precision  when  only  inertial  forcing  is  present. 
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Aerodynamic  Model 


As  discussed  in  the  literature  review,  a  number  of  researchers  have  investigated  the  aerody¬ 
namics  of  flapping  flight.  Most  have  used  quasi  steady  methods  [2,  7,  8,  26],  computational 
fluid  dynamics  (CFD)  [38,  10,  46],  or  a  hybrid  form  of  the  two.  The  purpose  of  this  thesis 
was  to  create  a  flight  dynamic  model  for  the  eventual  purpose  of  controller  synthesis  and  the 
analysis  of  flight  dynamic  metrics  such  as  stability  and  maneuverability.  The  aerodynamic 
models  discussed  above  are  inadequate  for  this  application.  The  quasi  steady  methods  often 
run  efficiently  but  are  not  able  to  fully  capture  the  complicated  phenomena  surrounding  the 
complex  flapping  motion.  CFD  methods  are  able  to  capture  these  phenomena  but  can  take 
days  to  run  a  single  period  of  analysis,  much  too  long  for  use  in  a  controller.  For  these 
reasons  it  was  important  to  create  a  reduced  order  model  capable  of  accurately  capturing 
these  phenomena  but  also  one  that  runs  efficiently  with  run-times  more  on  the  order  of  the 
quasi  steady  methods.  The  airloads  model  developed  by  Peters  et  al.  [15]  coupled  with 
an  assumed  inflow  model  allows  for  the  degree  of  accuracy  necessary  while  still  running  in 
reasonable  amounts  of  time  for  controller  synthesis. 
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3.1  Airloads  Model 

Peters’  airloads  model  allows  for  large  frame  motion  (such  as  the  prescribed  motion  of  the 
wings)  as  well  as  small  deformation  from  this  frame.  As  this  thesis  deals  only  with  rigid 
bodies  the  small  deformation  component  is  not  used.  The  airloads  model  also  provides  for 
inclusion  of  a  completely  separate  inflow  model  that  can  be  varied  independently  from  the 
airloads  model  itself.  This  allows  the  quasi  steady  airloads  model  to  include  as  detailed 
a  wake  (inflow)  model  as  desired,  up  to  and  including  a  fully  3D  unsteady  wake.  Thus 
by  carefully  calculating  and  including  all  effects  it  is  possible  to  take  what  begins  as  a  2D 
airloads  theory  and  create  a  fully  encompassing  3D  aerodynamic  model.  The  importance 
here  is  that  this  inflow  model  can  be  varied  to  be  as  accurate  and  as  efficient  as  desired,  all 
completely  separately  from  the  quasi  steady  airloads  model. 

The  airloads  model  itself  is  derived  from  potential  flow  and  begins  with  the  non-penetration 
boundary  condition  due  to  the  velocity  from  the  movement  of  the  airfoil,  the  freestream  flow, 
and  any  induced  flow  (inflow). 

dh  dh  x  .  . 

w  =  u°di+m+Vo  +  Vlb  (X1) 

where  w  is  the  total  induced  flow,  u0  and  v0  are  the  velocities  of  the  airfoil  in  the  x  and  y 
direction  (illustrated  in  Figure  3.1),  v\  is  the  velocity  gradient,  and  h  terms  represent  the 
small  deformation  from  the  large  frame  motion  which  in  the  present  work  are  set  to  zero. 

The  derivation  is  presented  with  more  details  in  Peters  et  al.  [15].  The  components  of 
the  total  induced  flow  are  related  to  the  general  airloads. 


T0  =  M0  (w0) 

(3.2) 

T\  =  bw0  +  UqWi 

(3.3) 

b 

T2  = 

(3.4) 

Ts  =  fT0 

(3.5) 
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Figure  3.1:  Peters’  airloads  model  coordinate  system 

where  /  is  the  reversed  flow  parameter.  These  equations  are  greatly  simplified  from  their 
original  form  due  to  the  assumption  that  all  bodies  are  rigid.  This  also  leads  to  wq  =  t’o  and 
W\  =  V\  (wm  =  0  where  m  A  2).  For  more  detailed  equations  with  deformation  terms,  see 
Peters  et  al.  [15].  The  velocities  of  the  airfoil  with  respect  to  the  air,  u0  and  v0  as  well  as  the 
velocity  gradient  V\,  are  written  in  the  wing  frame  and  are  dependent  on  the  fuselage  velocity, 
the  prescribed  wing  velocity  with  respect  to  the  fuselage,  any  freestream  flow  present,  and 
inflow  due  to  shed  vorticity. 

u0  =  -yk  cos (r)  +  yfi  sin(r)  -  ^A*  +  WVXB  (3.6) 

v0  =  —yfi  cos(r)  -  yksm(r)  -  w\y  +  WVZB  (3.7) 

vi  =  -hr  -  bwuby  (3.8) 

where  y  is  the  spanwise  location,  is  the  inflow  velocity  in  the  wing  x  direction,  w\y  is 
the  inflow  velocity  in  the  wing  z  direction,  UVB  is  the  velocity  of  the  body  in  the  wing  x 
direction,  WVB  is  the  velocity  of  the  body  in  the  wing  z  direction,  b  is  the  semi-chord,  and 
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wujy  is  the  angular  velocity  of  the  body  in  the  wing  y  direction.  The  general  airloads  can 
then  be  resolved  into  a  pressure  difference  A P  (again,  details  in  Peters  et  al.  [15])  which 
is  integrated  along  the  chord  to  get  the  aerodynamic  lift,  thrust,  and  moment  in  the  wing 
frame  (per  unit  length). 


L0  =  -2irpb  (/»(,'' Vi  +  i  (W'o  +  «o«i)j 


T  =  2npbfv% 

Li  =  77 pb  fuoVo  -  T  '1 1 


(3.9) 

(3,10) 

(3-11) 


These  forces  can  then  be  integrated  over  the  span  or  summed  using  blade  element  theory  in 
order  to  get  the  total  forcing  and  moment  for  the  entire  wing. 


3.1.1  Reversed  Flow 

Reversed  flow  is  the  term  used  to  refer  to  the  case  when  the  traditional  leading  edge  is  no 
longer  the  actual  leading  edge.  In  other  words,  what  is  usually  the  trailing  edge  is  leading 
with  respect  to  the  motion  of  the  airfoil.  Several  studies  have  been  done  that  illustrate  the 
effectiveness  in  using  reversed  flow  (through  timing  of  wing  rotation  with  respect  to  stroke) 
as  a  control  mechanism  [2,  30,  14].  The  rotation  of  the  wing  can  be  altered  to  change  the 
reaction  forces  caused  by  the  rotational  effects.  If  the  wing  flips  early  (before  the  end  of 
the  half  stroke)  then  the  resulting  total  force  should  be  upward.  If  the  flip  is  late  then  the 
resulting  force  (due  to  rotational  effects)  would  be  a  downward  force.  If  the  flip  spans  both 
half  strokes  then  the  resulting  downward  and  upward  forces  (due  to  rotational  effects)  cancel. 
In  Figure  3.2  you  can  see  how  the  advanced  or  delayed  rotation  can  add  or  subtract  lift  from 
the  case  where  rotation  is  symmetric.  It  is  also  clear  that  the  trailing  edge  leads  at  some 
points  in  the  stroke.  Therefore  proper  timing  of  wing  rotation  on  each  wing  can  create  a 
force  imbalance  that  can  be  used  for  control  [2],  The  reversed  flow  caused  by  these  varying 
rotations  (anytime  the  trailing  edge  is  leading)  is  accounted  for  the  the  variable  /  in  Equation 
3.5.  This  parameter,  /,  is  called  the  reversed  flow  parameter  and  ensures  that  the  A P  =  0 
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advanced 

■M —  downstroke  s 


symmetrical 

—  downstroke 


\ 
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Figure  3.2:  Results  of  alteration  to  wing  flip  timing  [2] 


at  the  true  trailing  edge.  There  are  several  cases  for  /: 

/  =  1  (reversed  flow  neglected)  (3-12) 

/  =  —  1  (trailing  edge  always  at  what  is  normally  the  leading  edge)  (3.13) 

/  =  T"~~T  (full  reversed  flow)  (3.14) 

/  =  cos(a)  (soft  reversed  flow)  (3.15) 


where  a  is  the  airfoil  angle  of  attack  as  defined  by  Peters  [15].  In  the  case  of  full  reversed 
flow  the  sign  of  /  changes  with  the  sign  of  uq  moving  the  stagnation  point  and  thus  the  lift 
direction.  Soft  reversed  flow  makes  a  smooth  transition  between  these  two  cases  and  makes 
lift  proportional  to  blll(22a'>  instead  of  sin(a). 


3.2  Assumed  Inflow  Models 

The  quasi-steady  airloads  model  can  be  coupled  with  an  inflow  model  to  account  for  the 
unsteady  effects  present  due  to  the  complicated  nature  of  flapping  flight.  These  inflow  mod- 
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els  are  completely  separate  from  the  airloads  model  as  shown  by  the  separate  A  terms  in 
Equations  3.6  and  3.7.  Thus  different  inflow  models  can  be  developed  and  applied  to  the 
aerodynamic  forcing  equations  allowing  for  a  wide  range  of  accuracy  and  efficiency  (run 
time).  Theoretically,  given  a  fully  3D  inflow  that  takes  into  account  all  the  wake  effects  at 
all  times,  one  could  create  a  perfect  3D  aerodynamic  code  given  the  inflow  model  and  the 
2D  airloads  theory.  The  inflow  models  discussed  below  were  used  in  various  stages  of  this 
thesis  work. 


3.2.1  Momentum  Disk  Theory 


Helicopter  theory  often  makes  use  of  inflow  in  calculating  aerodynamic  loads  as  well  as 
calculating  power  and  efficiency.  A  common  initial  approach  to  inflow  when  dealing  with 
helicopters  is  to  assume  a  constant  inflow  based  on  momentum  disk  theory  as  detailed  by 
Johnson  [47].  Using  momentum  theory  the  thrust  generated  by  a  helicopter  rotor  can  be 
related  to  the  induced  velocity  by  the  equation: 


vh  = 


(3,16) 


where  is  the  inflow,  T  is  the  thrust  from  the  rotor,  p  is  the  density  of  the  air,  and  A  is  the 
area  swept  by  the  rotor.  The  thrust  necessary  to  counteract  gravity  at  trim  is  just  mtotaig- 
The  area  swept  by  the  wings  is  2nmaxR2LU.  This  leads  to  the  inflow  at  hover  of: 


A  = 


TFltotalQ 

^P^rnaxFiw 


(3.17) 


Alternatively  one  can  calculate  the  lift  force  in  the  body  frame  without  including  inflow  and 
use  the  average  lift  over  a  period  as  the  thrust  from  the  rotor,  T. 


A  =\l  ybFA)*  (3.18) 

V  4 pKmaxRw 

Once  a  new  lift  force  is  calculated  using  this  inflow  it  can  be  substituted  as  a  new  rotor 
thrust  force  and  the  process  iterated  to  a  specified  convergence  condition. 
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In  order  to  be  used  in  Equations  3.6  and  3.7  this  velocity  must  be  rotated  into  the  wing 
frame  and  the  correct  components  used. 

{A,,  Xy,  Az}  T  =  [Tws]s  {0, 0,  A}  T  (3.19) 


3.3  Power  Calculations 

The  power  is  calculated  by  using  Equations  3.9,  3.10,  and  3.11  for  each  blade  element  (where  r/ 
is  the  blade  width)  and  multiplying  them  by  the  wing  velocity  in  the  corresponding  direction. 


Power  =  Tuo  +  L0vq  +  L\V\ 


(3.20) 


Expanding  this  equation  and  simplifying  terms  yields  the  equation  below: 


Power  =  npbrj 


-bi! 0v0  -  -ViVi 


(3.21! 


Because  the  velocity  functions  are  periodic  they  can  be  expanded  using  Fourier  series 
(using  vq  without  inflow  as  an  example): 


OO  OO 

VqVq  =  EE  mu  [an  cos(nut )  +  bn  sin(naA)]  [— am  sin (mut)  +  bm  cos (mut)]  (3.22) 

ri—l  m=  1 

Integrating  these  equations  over  one  period  to  get  the  average  power  leads  to  cancellation 
of  the  terms.  When  n  is  not  equal  to  m  all  terms  cancel  as  they  are  orthogonal.  When  n  is 
equal  to  m: 


cos  (nut)  sin  {mut)dt  = 


'  period 


/  period 


cos  (nut)  sin  {nut)dt  =  0  (orthogonal)  (3.23) 
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anbm  /  cos(noA)  cos(mut)dt  —  ambn  /  sm(nut)  sin(mut)dt  (3.24) 

J  period  J  period 


anbn  /  cos2  {nut)dt  —  anbn  /  sin 2(ncut)dt 

J  period  J  period 


(3.25) 


=  anbn  /  [cos2  (nut)  —  sin2(naA)]  dt 

J  period 


(3.26) 


n  n  b7l 


f 

1  +  cos  (2  nwt) 

1  —  cos(2na;t ) 

J  period 

2 

2 

dt 


(3.27) 


=  anbn  /  cos(2 nut)dt  (which  is  zero  when  the  period  is  equal  to  — )  (3.28) 

J period  ^ 


The  solution  for  i’i  is  of  the  same  form,  showing  that  for  periodic  velocities  without  inflow 
the  average  power  is  always  zero.  This  is  significant  in  that  if  one  were  to  optimize  for  average 
power  there  would  not  be  a  solution,  emphasizing  how  important  the  inflow  solution  is  to 
this  problem.  An  accurate  inflow  solution  is  essential  when  considering  matters  of  power.  If 
the  inflow  is  included  in  the  power  calculations  the  resulting  power  is  equal  to  that  shown  in 
Equation  3.29. 

Power  =  npbr]2fu0v0Xy  +  2/u0A2  -  2  fv%\x  -  2fv0XxXy  +  u0v1Xy  +  XxXyvi  -  \vxVi  (3.29) 


3.4  Validation 

For  verifying  the  aerodynamic  model  a  comparison  was  made  with  the  quasi  steady  aero¬ 
dynamics  used  in  Stanford  [3].  Stanford  [3]  used  the  2D  quasi-steady  model  developed  by 
Berman  and  Wang  [26].  The  aerodynamics  of  Berman  and  Wang  [26]  are  applicable  only 
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to  rigid  airfoils  undergoing  large  motion.  Peters  aerodynamics  can  handle  small  flexible  mo¬ 
tions  in  addition  to  the  large  rigid  motion  but  the  flexible  part  is  not  considered  in  this  work. 
There  are  several  key  differences,  however,  that  must  be  accounted  for  before  the  models  can 
be  directly  compared.  Examining  the  equations  used  for  the  forcing  in  Berman  and  Wang 
[26]: 


Fy'e 


fv  = 


m22  ■  vz>e  ■  ip  +  pf  ■  T  •  Vgf  -  mn  •  ay 


-mu  •  Vy>e  ■  Ip  +  Pf  •  r  •  Vy'e  -  m22  •  az>e 


T  =  —0.5  •  Ct  ■  c  ■  \  v2,  +  v2,  ■  sin(2  •  a)  +  0.5  •  Cr  ■  c2  ■  ip 

*'  Ue 


dFZ  1 

ye 

dx'e  _ 

■  dx'e 

(3.30) 

dFvz, 

ze 

dx 

■  dx 

(3.31) 

Cr  ■  c 

2 

(3.32) 

dx ' 


\  Fv.  1 

|  =  0.5-pf-c-(CD(0)  ■  cos2 (a)  +  CD( 7t/2)  •  sin2(o:))  ■  yjv2,  +  v2z,  -  j 

I  vv'e  \ 

l  J 

{  >’4  | 

(3.33) 


Which  when  converted  to  the  variables  used  in  this  work  yield: 


Ln  — 


1C 


-2 CTpb  (  u0v0f  +  ^-£T  (^0  +  u0v i)  j  r]  ^ 


dF» 

Asp 


(3.34) 


dFv, 

T  =  2CTpbfvlr)  +  CRpbv0Vir]  - 


(3.35) 


and  f  =  — .  °  (3.36) 

Equations  3.34,  3.35,  and  3.36  illustrate  the  key  differences  in  the  aerodynamics  being 
used  for  this  work  and  those  presented  by  Berman  [26].  The  last  terms  in  Equations  3.34  and 
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Figure  3.3:  Comparison  between  altered  Peters’  airloads  model  and  Stanford  [3]  results 


3.35  represent  viscous  drag  terms  and  for  comparison  are  calculated  using  the  same  method 
presented  in  Equation  3.33.  The  Ct  and  Cr  presented  in  these  equations  are  empirically 
determined  coefficients  that  can  be  set  to  n  to  agree  with  Peters  theory  or  the  coefficients 
from  Berman  [26]  for  comparison  to  that  work.  As  seen  in  Equation  3.36  the  reversed  flow 
coefficient  must  be  set  to  the  soft  reversed  flow  condition  in  order  to  match  the  blll^'Q^  trend 
from  Equation  3.32.  The  study  also  uses  an  elliptic  wing  as  defined  in  Equation  3.37  where  c 
is  the  chord,  c  is  the  average  chord,  y  is  the  span  location,  and  Rw  is  the  overall  wing  length. 
The  last  major  difference  is  the  presence  of  a  VqV\  term  in  the  thrust  where  there  is  not  one  in 
Peters  theory.  This  term  turns  the  force  due  to  the  u0v i  term,  usually  perpendicular  to  the 
wing,  and  makes  it  perpendicular  to  the  flow.  Also  ignored  here  are  the  added  mass  terms  in 
the  chord-wise  direction.  When  these  terms  are  set  to  agree  (viscous  drag  added,  coefficients 
of  thrust  and  rotation  equal,  added  mass  ignored,  soft  reversed  flow  added,  elliptic  wing,  and 
VqVi  term  added)  the  results  are  as  shown  in  Figure  3.3  and  exhibit  exact  agreement  (the 
effect  of  each  of  the  various  assumptions  is  considered  further  in  Section  5.2).  The  parameters 
used  in  this  comparison  are  shown  in  Table  3.1. 

(3.37) 


4c 


c  = 


7 r 


i  _  jr. 
1  R? 

j. 
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Parameter 

Value 

Units 

V 

26 

Hz 

Ryj 

50  ■  10"3 

m 

K0 

0 

radians 

K'max 

7t/3 

radians 

A) 

tt/2 

radians 

7 max 

7t/4 

radians 

&tau 

0 

radians 

IM) 

0 

radians 

l^max 

0 

radians 

<t>R 

0 

radians 

Ct 

1.833 

- 

Cr 

7 r 

- 

cD{  0) 

0.21 

- 

CD(  tt/2) 

3.35 

- 

c 

18.2  •  10“3 

m 

rriR 

1.648  ■  10“3 

kg 

mw 

4.7-  10“5 

kg 

Table  3.1:  Parameters  used  in  aerodynamic  validation  (Figure  3.3) 
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Trim  and  Stability 

4.1  Stroke- Averaged  System 

Since  the  wing  frequencies  are  much  higher  than  the  flight  dynamic  frequencies,  flapping  wing 
aerodynamics  on  a  MAV  scale  are  often  approximated  on  a  stroke-averaged  basis  as  these 
are  sufficiently  accurate  enough  on  the  time  scales  of  the  rigid-body  dynamics  [34,  36].  Here 
the  linearized  forcing  is  determined  using  perturbations  from  a  reference  state  to  determine 
•Jacobians  and  initial  forcing  vectors.  In  this  particular  study  only  symmetric  variables  are 
considered. 

The  system  equations  are  of  the  form  (Subsection  2.2.2): 

[M(((t))]x  =  f(x,((t))  (4.1) 

One  can  average  the  forcing  over  one  flapping  cycle  for  a  specified  periodic  kinematics 
C(t)  and  constant  prescribed  state  x0.  This  yields  an  average  forcing  dependent  on  kinematic 
variables,  Co,  and  the  state.  The  kinematic  variables,  Co  are  variables  that  parametrize  the 
actual  kinematics,  e.g.,  kq,  Kmax,  rmax,  0r. 
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(4.2) 


It  is  also  necessary  to  average  each  entry  of  the  mass  matrix,  M(£(t)),  for  one  stroke 
period.  This  leads  to  the  stroke  averaged  system: 


favgip^i  Co) 


(4.3) 


4.2  Stroke-Averaged  Nonlinear  Trim 

The  objective  of  this  section  is  to  determine  the  trim  conditions  of  the  stroke-averaged  system 


by  equating  the  stroke-averaged  forces  to  zero.  This  leads  to  no  change  in  the  velocity  states 
of  the  system  over  a  period.  Given  the  stroke-averaged  forcing  (Equation  4.3)  it  is  possible 
to  solve  for  a  stroke  averaged  non-linear  trim  condition.  In  this  case,  because  the  system 
has  been  stroke-averaged,  the  goal  of  the  trim  analysis  is  to  get  x  =  0  which  is  equivalent 


to  favg  —  0.  Since  the  aerodynamic  forces  are  a  nonlinear  function  of  the  state  variables 
as  well  as  the  control  variables,  a  nonlinear  iterative  solution  technique  is  required.  Here  a 


Newton-Raphson  method  is  used.  The  Jacobian  (derivative  of  the  function  with  respect  to 


the  variables)  is  calculated  using  finite  difference. 

The  first  set  of  trim  variables,  £,  considered  are  those  of  the  symmetric  states  relevant  to 
the  aerodynamics  and  gravity,  u,  w,  q,  and  0  (while  holding  the  wing  kinematics  constant). 
First  it  is  necessary  to  calculate  the  initial  stroke-averaged  aerodynamic  and  gravitational 
forcing  vectors,  Fa0  and  FGo  (consisting  of  the  forward  forcing  X,  the  vertical  forcing  Z,  and 
the  pitch  moment  M,  as  well  as  a  zero  in  the  fourth  row  from  the  fourth  equation  q  =  6  =  0) 
based  on  the  trim  variables.  Perturbing  each  of  the  entries  in  the  vector  £  and  comparing  the 
period  averaged  forcing  to  the  initial  values,  Fa0  and  Fg  o,  provides  the  the  columns  of  the 
Jacobians,  ^Fa.  anci  ^Fa.  These  Jacobians  are  in  effect  the  stability  derivatives  and  would 
correspond,  for  example,  to  the  slopes  of  the  curves  found  in  Figures  5.6,  5.7,  5.8,  and  5.9  if 
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the  initial  state  were  chosen  to  be  the  zero  state  for  each  of  the  variables.  This  provides  the 
average  forcing  linearized  about  the  trim  variables: 


favg(x,  Co) 


ofa 

dx 


(AC}  + 


ofg 

dx 


{Aft  +  {FAg}  +  {FGo} 


(4.4) 


It  is  then  possible  to  solve  for  a  new  initial  state  by  iterating  the  equation: 


Ci  —  Co 


8Fa  dFr 


-\  -i 


(Fa  o  +  FGo) 


(4.5) 


L  d£  d£ 

Iteration  continues  until  the  error  metric  (error  vector  squared)  is  less  than  1CT10.  In 
order  to  prescribe  certain  motions  such  as  hover,  forward  flight,  and  climb,  more  constraints 
are  needed  in  the  trim  analysis.  In  particular,  the  derivatives  of  the  vertical  and  horizontal 
inertial  positions  have  to  be  prescribed  ( iXs  and  jig).  This  adds  two  more  equations  to  the 
system,  namely  Equations  4.6  and  4.7. 


iXb  =  prescribed 

(4.6) 

iZs  =  prescribed 

(4.7) 

This  necessitates  the  use  of  all  four  symmetric  states  and  two  kinematic  variables  as  trim 
variables.  The  kinematic  variables  chosen  for  this  study  as  control  variables  are  the  stroke 
offset,  k0,  and  the  maximum  magnitude  of  rotation,  Tmax.  Otherwise  the  trim  solution  is  the 
same  as  given  in  Equations  4.4  and  4.5. 


4.3  Stroke- Averaged  Linearized  Stability  About  Trim 

After  a  trim  state  has  been  calculated  it  is  possible  to  calculate  the  stability  about  this 
trim  state  by  simply  taking  the  eigenvalues  of  the  linearized  stability  matrix,  +  ^r\ , 
calculated  about  the  trim  state.  The  above  matrix  is  the  subset  of  the  matrix  computed  for 
the  trim  solution 
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4.4  Periodic  Trim 

Without  stroke  averaging,  the  system  (Equation  2.43)  is  a  time-varying  system.  For  periodic 
kinematics,  the  system  becomes  a  time- varying  periodic  system.  A  periodic  system  can  have 
a  periodic  trim  solution  which  can  calculated  using  a  periodic  trim  solution  technique.  The 
goal  of  a  periodic  shooting  calculation  is  to  find  the  initial  condition  which  leads  to  a  periodic 
solution.  The  simulation  of  the  system  over  one  period  using  the  initial  condition  results  in 
a  final  state  equal  to  the  initial  one. 


x(t0,  x0)  =  x(t0  +  T,  x0)  (4.8) 

As  this  analysis  focuses  on  the  symmetric  variables  only  the  longitudinal  states  are  of 
concern.  Because  of  this  only  two  kinematic  variables  are  needed  for  control,  as  stated  above. 
In  order  to  solve  for  the  initial  trim  state  a  periodic  shooting  method  is  implemented  whereby 
each  of  the  symmetric  states  is  perturbed  and  time  marched  through  one  period  of  motion 
and  then  compared  to  the  unperturbed  final  state.  Stated  explicitly,  after  time  marching  the 
initial  guess  through  one  period  the  difference  in  the  final  and  initial  states  is  the  error: 

x9t  —  x9q  =  error  (4.9) 

Where  xaT  is  x(T,  x90)  (4.10) 

The  error  vector  in  this  case  is  a  non-linear  function  of  the  trim  variables,  £.  Squaring 
the  error  vector  gives  the  error  metric  by  which  convergence  is  determined. 

Starting  with  an  initial  guess  for  the  trim  variables  £o,  one  can  iterate  using  the  Newton 
Raphson  scheme  to  get  a  trim  condition  (when  the  error  metric  is  less  than  10~18): 

I0  +  |2A?_t  +  |iA5  (4.11) 
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6  =  £o  + 


(4.12) 


It  is  important  to  note  that,  except  for  the  case  of  completely  prescribed  kinematics,  extra 
equations  (beyond  the  four  state  equations)  are  needed  to  prescribe  the  motion  of  the  system 
as  desired.  In  order  to  prescribe  hover,  forward  flight,  and  climb,  the  averages  of  the  deriva¬ 


tives  of  the  vertical  and  horizontal  inertial  positions  have  to  be  prescribed  (iX b avg  and i^B avg)  ■ 
This  adds  two  more  equations  to  the  system,  namely  Equations  4.13  and  4.14. 


(4,13) 


(4,14) 


This  necessitates  the  use  of  all  four  symmetric  states  and  two  kinematic  variables  as  trim 
variables  (same  trim  variables  as  the  stroke-averaged  case). 

4.5  Stability  of  the  Periodic  System  (Floquet  Analysis) 

To  calculate  the  stability  of  the  system  we  calculate  the  state  transition  matrix.  The  state 
transition  matrix  gives  the  effect  of  perturbation  (about  the  trim)  in  the  initial  state  on  the 
state  at  the  end  of  one  period.  If  the  perturbation  grow  then  we  have  instability  otherwise 
we  have  a  stable  system.  In  order  to  do  this  we  perturb  each  of  the  states  in  the  initial 
condition,  xo,  and  use  the  perturbation  in  the  resulting  solutions  at  the  end  of  one  period 
to  form  the  columns  of  the  state  transition  matrix  0.  Note  that  this  matrix  is  a  subset  of 
the  matrix  calculated  for  trim  calculations  based  on  periodic  shooting.  For  stability  we  have 
0  =  eAT  where  the  eigenvalues  of  A  determine  the  stability  of  the  system. 
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5.1  System  Parameters 

The  system  parameters  used  in  this  study  are  altered  Manduca  sexta  parameters  derived 
from  [48,  49,  50]  and  are  listed  below  in  Table  5.1.  The  vehicle  itself  is  not  meant  to  be  an 
exact  representation  of  the  Manduca  sexta  but  only  based  on  one,  and  so  the  values  may  be 
slightly  different  than  those  found  in  literature.  There  are  also  crucial  differences  in  that  the 
body  is  modeled  as  an  ellipsoid  and  the  wings  are  rectangular.  These  parameters  may  also 
be  varied  in  specific  portions  of  this  work  as  described  throughout. 
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Symbol 

Parameter 

Value 

Units 

mB 

Body  mass 

1.780. 10"3 

kg 

mw 

Wing  mass 

2.133. 10'5 

kg 

Tb 

1B 

Body  inertia  tensor 

Equation  5.1 

kg  .  m2 

tW 

1W 

Wing  inertia  tensor 

Equations  5.2,  5.3 

kg  .  m2 

rBH 

Vector  from  the  body  CG  to  the  hinge  point 

Equations  5.4,  5.5 

m 

rHW 

Vector  from  the  hinge  point  to  the  wing  CG 

Equations  5.6,  5.7 

m 

C 

Average  chord 

0.019427 

m 

Rw 

Wing  length 

0.05055 

m 

sw 

Wing  area 

0.98203 . 10"3 

2 

m 

V 

Flapping  frequency 

25 

Hz 

Table  5.1:  Morphologic  parameters 


Assuming  an  ellipsoid  body  with  body  length  (46.377  mm)  and  radius  (5.93  mm)  from 
[48,  49]  provides  a  body  inertia  tensor  as  shown  in  Equation  5.1. 


tb  - 

b1B  ~ 


10  8  kg  .  m2 


(5.1) 


77.821  0  0 

0  77.821  0 

0  0  2.5307 

Assuming  rectangular  wings  using  the  parameters  shown  in  Table  5.1  provides  the  right 
and  left  wing  inertia  tensors. 

4.5421  0  0 

0  0.67084  0 

0  0  5.2129 


TWr  _ 

Wr1Wr 


' 1  9  kg  . m2 


(5.2) 


TWi  _ 
wi1Wi  ~ 


4.5421  0  0 

0  0.67084  0 

0  0  5.2129 


10  9  kg  .  mz 


(5.3) 
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No  data  was  found  for  the  vectors  from  the  body  CG  to  the  hinge 

Equations  5.4  and  5.5  were  assumed  where  appropriate. 

so  the  values  in 

brBHr  =  [0j  5.58,  — 10]T  .  10"3m 

(5,4) 

brBHl  =  [0,  -5.58,  -10]r  .  10“3m 

(5,5) 

Assuming  rotation  about  the  midchord  and  the  rectangular  wing  profile  provides  the 

following  hinge  to  wing  CG  vectors  for  the  right  and  left  wings. 

wrHWr  =  [0,  0.025275,  0]Tm 

(5.6) 

wrHWr  =  [0,  -0.025275,  0]T  .  10~3rn 

(5.7) 

5.2  Effects  of  Model  Assumptions  on  Aerodynamic  Forc¬ 
ing 

As  mentioned  in  Sections  3.2  and  3.4  there  are  many  different  modeling  considerations  when 
trying  to  accurately  model  the  complex  flapping  wing  system  of  a  flapping  wing  MAV.  In 
this  section  the  effect  of  some  of  these  models  on  the  aerodynamic  forcing  will  be  examined. 
The  parameters  used  were  the  system  parameters  described  above  with  kinematics  as  shown 
in  Table  5.2. 
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Parameter 

Value  (degrees) 

K  0 

0 

fc'max 

0.9594 

To 

vr/2 

7 max 

1.4667 

(j)T 

0 

fj-o 

0 

H'max 

0 

AU 

0 

Table  5.2:  Kinematic  Parameters 


5.2.1  Effect  of  Reversed  Flow 


Figure  5.1:  Comparison  of  aerodynamic  forcing  curves  showing  the  effect  of  adding  reversed  flow 

Figure  5.1  shows  the  effect  of  adding  soft  reversed  flow  (defined  in  Section  3.1.1)  to  the 
quasi-steady  aerodynamic  model.  Examining  the  body  frame  forcing  curves  in  Figure  5.1 
shows  that  adding  reversed  flow  does  not  affect  the  body  frame  thrust  force  in  an  average 
or  an  instantaneous  sense  but  it  does  reduce  the  magnitude  of  the  body  lift  force  at  all 
times  during  the  period.  This  is  to  be  expected  given  that  the  inclusion  of  soft  reversed  flow 
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switches  the  trend  of  lift  from  sin(a)  to  The  soft  reversed  flow  also  explains  the  double 

peaks  per  half  stroke  in  the  body  lift  curves.  As  the  wing  rotates  down  from  ninety  degrees 
the  lift  increases  up  until  the  point  it  reaches  forty-five  degrees  at  which  point  it  begins  to 
reduce  again,  at  the  middle  of  the  half  stroke  it  reaches  a  minima  before  increasing  its  as 
it  increases  in  angle  of  attack  again.  Towards  the  end  of  the  half  stroke  the  lift  magnitude 
again  reduces  as  the  angle  of  attack  moves  from  forty-five  degrees  back  to  ninety.  This  same 
pattern  repeats  on  the  other  half  stroke. 

5.2.2  Effect  of  Empirical  Coefficients 


Figure  5.2:  Comparison  of  aerodynamic  forcing  with  the  addition  of  empirical  constants 

Figure  5.2  shows  the  results  of  adding  the  empirical  coefficients  CT  =  1.678  shown  in  Equa¬ 
tions  3.34  and  3.35  to  the  results  obtained  in  Section  5.2.1.  Adding  the  empirical  coefficient 
into  the  forcing  model  again  does  not  have  an  effect  on  the  body  frame  thrust  force  but  does 
have  a  significant  affect  on  the  body  lift  force.  Decreasing  the  coefficient  of  thrust  relative  to 
the  coefficient  of  rotation  puts  more  emphasis  on  the  rotational  portions  of  the  wing  frame 
lift  and  thrust  equations  and  decreases  the  overall  magnitude  of  the  forcing. 
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5.2.3  Effect  of  Viscosity 


Figure  5.3:  Comparison  of  aerodynamic  forcing  with  the  addition  of  viscosity 

Figure  5.3  shows  that  effect  of  adding  the  viscous  forcing  shown  in  Equations  3.34  and  3.35  to 
the  results  obtained  in  Section  5.2.2.  The  addition  of  viscosity  changes  only  the  instantaneous 
body  thrust.  It  increases  drag  for  both  the  forward  and  backward  stroke  increasing  the  power 
required.  It  does  not  affect  the  stroke-averaged  thrust  as  its  effects  in  the  body  frame  due  to 
the  forward  stroke  are  canceled  by  those  on  the  back  stroke. 

5.2.4  Effect  of  vqV\  Term 


Figure  5.4:  Comparison  of  aerodynamic  forcing  curves  with  the  addition  of  the  vqV\  term 
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Figure  5.4  shows  the  effect  of  adding  the  v0vi  term  shown  in  Equation  3.35  to  the  results 
obtained  in  Section  5.2.3.  It  shows  that  while  adding  the  vov±  term  to  Equation  3.35  does 
affect  instantaneous  forcing  values  in  the  body  frame  it  does  not  affect  the  overall  averages 
over  the  period. 

5.2.5  Inflow 


Figure  5.5:  Comparison  of  aerodynamic  forcing  curves  with  the  addition  of  inflow 

Figure  5.5  shows  the  effect  of  adding  the  prescribed  inflow  as  discussed  in  Section  3.2  to  the 
results  obtained  in  Section  5.2.4.  It  can  be  seen  that  adding  a  constant  inflow  model  does 
not  affect  the  overall  average  of  the  body  thrust  but  does  affect  both  body  thrust  and  lift 
instantaneously  and  reduces  the  stroke  average  of  the  body  lift  force  due  to  the  energy  lost 
in  the  wake. 


5.3  Nominal  System  Dynamics 

The  nominal  case  consists  of  a  system  utilizing  the  system  parameters  discussed  in  Section 
5.1.  The  aerodynamics  utilize  soft  reversed  flow  but  do  not  include  viscous  drag  or  inflow. 
All  time  marching  was  accomplished  using  the  MATLAB  program  ode45  using  a  relative 
tolerance  of  1CT10  and  an  absolute  tolerance  of  1CT15. 
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5.3.1  Influence  of  State  Variables  on  Forcing 

The  effect  of  state  and  potential  kinematic  control  variables  (particularly  k0,  nmax ,  To?  and 
Tmax)  is  examined  over  the  course  of  a  stroke  cycle  (stroke-averaged  forcing).  The  slopes  of 
the  resulting  data  can  be  interpreted  as  aerodynamic  (stability  and  control)  derivatives.  Both 
symmetric  and  antisymmetric  variations  are  considered  though  it  is  important  to  note  that 
the  kinematics  and  geometry  are  symmetric  so  changes  in  the  symmetric  state  variables  can 
only  lead  to  symmetric  load  variations  (though  not  necessarily  vice  versa).  Using  the  nominal 
case  (reversed  flow,  no  empirical  coefficients,  no  inflow,  no  viscous  drag,  no  deviation,  AC 
above  the  body  CG),  and  a  hover  trim  condition,  the  results  are  as  shown  below  in  Figures 
5.6,  5.7,  5.8,  and  5.9.  For  display  purposes  the  velocities,  forces,  and  moments  are  normalized 
by  the  characteristic  velocity,  force,  and  moment  as  described  in  Equations  5.8,  5.9,  and  5.10. 


1  2Rwu  sin(ftmax) 


(5.8) 


F  =  pV2Sw 


(5.9) 


M  =  pV2Svlc 


(5,10) 


Figure  5.6:  Influence  of  longitudinal  variables  on  longitudinal  forcing 


Figure  5.6  shows  the  variation  in  the  symmetric  aerodynamic  loads  (x  direction  force 
XB ,  z  direction  force  ZB ,  and  pitching  moment  MB )  due  to  perturbation  in  the  symmetric 
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velocities.  The  velocity  perturbation  in  the  x  direction,  uB ,  leads  to  little  change  in  loads 
other  than  the  pitching  moment,  MB .  Due  to  the  lack  of  any  viscous  drag  forces  there  is  no 
change  in  x  direction  force  (as  seen  later),  force  in  the  z  direction  experiences  a  similar  lack 
of  change  because  increase  in  the  force  from  increased  velocity  in  one  direction  are  offset  by 
decrease  in  forcing  while  flapping  in  the  opposite  direction.  The  significant  moment  effect  is 
due  to  the  fashion  in  which  2-D  pitching  moment,  Li,  is  calculated  in  Equation  3.11.  The 
uqVo  term  leads  to  an  additive  effect  on  the  forward  and  backward  stroke,  for  example,  if  the 
perturbation  is  positive  then  the  increase  in  these  two  quantities  on  the  forward  stroke  leads 
to  an  increased  pitching  moment.  On  the  reverse  stroke  a  decrease  in  both  of  these  terms 
means  that  the  L\  quantity  does  not  cancel  and  a  net  pitch  up  results  (after  rotation  into 
the  body  frame). 

When  perturbing  velocity  in  the  z  direction  the  only  significant  affect  is  a  change  in  the 
vertical  force  ZB .  An  increase  in  upward  velocity  lowers  the  effective  angle  of  attack  and 
thereby  decreases  the  vertical  force,  effectively  acting  as  a  damper  for  the  motion. 

Perturbing  the  pitch  rate  produces  a  significant  pitch  damping  as  well  as  a  change  in 
the  x  direction  forcing.  Perturbing  pitch  in  a  positive  direction  increases  the  angle  of  attack 
when  the  wings  are  rear  of  the  CG  and  decreases  it  when  forward  of  the  CG.  This  effectively 
produces  a  pitch  down  moment,  damping  out  the  pitch  perturbation  (the  opposite  of  this 
effect  is  present  with  a  negative  pitch  perturbation,  again  damping  out  the  motion).  The  x 
direction  forcing  perturbation  arises  from  the  uqV\  term  found  in  Equation  3.3.  Without  an 
equivalent  vovi  term  as  found  in  Berman  and  Wang  [26]  the  forcing  vector  that  results  from 
this  UqV i  term  is  perpendicular  to  the  wing  rather  than  to  the  flow.  Therefore  effective  angle 
of  attack  changes  result  in  forward  forcing  when  pitch  is  perturbed. 
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r  perturbation  x  10-4 


Figure  5.7:  Influence  of  lateral  variables  on  longitudinal  forcing 


Figure  5.7  shows  the  effect  of  lateral  motion  on  longitudinal  forces  ( y  direction  force  Y, 
roll  moment  L,  and  yaw  moment  N ),  especially  the  z  direction  force.  These  effects,  however, 
appear  to  be  small  when  compared  with  those  shown  in  Figures  5.6  and  5.8.  Particularly 
the  changes  in  forces  and  moments  due  to  roll  are  very  small  and  could  be  neglected. 


Figure  5.8:  Influence  of  lateral  variables  on  lateral  forcing 


Figure  5.8  shows  the  effect  of  lateral  velocity  perturbations  on  the  lateral  aerodynamic 
loads.  A  positive  perturbation  in  the  y  velocity  leads  to  a  significant  negative  rolling  moment. 
This  moment  is  due  to  the  increase  in  effective  angle  of  attack  on  the  right  wing  (the  wing 
toward  which  the  body  is  moving)  and  the  corresponding  decrease  in  angle  of  attack  of  the 
left  wing.  This  leads  to  a  net  negative  rolling  moment.  The  lack  of  a  yaw  as  well  as  side 
force  damping  is  due  to  the  absence  of  viscous  drag  force  in  this  model. 

A  significant  roll  damping  is  present  while  all  side  force  and  yaw  damping  is  negligible. 
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Figure  5.9:  Influence  of  control  variables  on  longitudinal  forcing 


In  Figure  5.9  the  effect  of  the  two  control  variables  demonstrates  how  the  vehicle  can  be 
controlled  for  symmetric  motion,  kq  perturbations  yield  pitching  moments  because  the  mean 
location  of  forcing  is  shifted  forward  or  back  from  the  CG  (positive  k0  shifts  the  wings  center 
of  motion  further  back).  This  allows  for  pitch  control,  which  also  yields  inertial  x  direction 
control  as  the  vehicle  can  be  tipped  forward  or  backward  and  then  z  body  direction  forces 
increased  or  decreased  to  achieve  inertial  x  direction  movement.  The  change  in  z  direction 
forcing  is  achieved  by  altering  the  rmax  kinematic  variable.  This  changes  the  effective  angle 
of  attack  of  the  wing. 


5.3.2  Stroke- Averaged  Nonlinear  Trim 

The  trim  condition  resulting  from  using  prescribed  kinematics  (k0  =  0,  Kmax  =  0.95944,  r0  = 
7t/2,  rmax  =  7t/4,  and  no  deviation)  are  shown  in  Table  5.3.  As  seen  in  the  table,  there  are 
at  least  four  possible  trim  conditions  which  where  found  by  starting  the  iterative  nonlinear 
solution  at  different  initial  guesses. 
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State 

Trim  Value 

u 

5.9715e-17 

w 

-7.2844 

q 

0 

0 

3.1416 

State 

Trim  Value 

u 

-1.7759e-15 

w 

-3.6136 

q 

0 

0 

1.4234e-15 

State 

Trim  Value 

u 

-2.4840e-16 

w 

1.7177 

q 

0 

0 

-1.2207e-16 

State 

Trim  Value 

u 

2.4680e-15 

w 

4.7524 

q 

0 

0 

3.1416 

(b)  (c)  (d) 


Table  5.3:  Multiple  sets  of  trim  variables  for  prescribed  kinematics 


In  Table  5.3  it  is  clear  from  the  first  set  of  data  that  in  order  to  maintain  a  constant 
state  with  the  given  kinematics,  the  vehicle  enters  a  constant  climb  in  order  to  counteract 
the  extra  lift  generated  by  the  wings.  The  other  states  remain  at  basically  zero  with  no 
forward  motion  or  pitching.  In  the  third  set  of  data  the  vehicle  is  flipped  over  as  seen  from 
the  0  value.  This  state  corresponds  to  a  flipped  over  constant  descent  which,  while  not  that 
interesting  physically,  illustrates  a  second  trim  condition  for  the  same  kinematic  variables. 
Two  more  possible  trim  states  are  shown  in  set  (b)  and  (d).  The  multiple  trim  conditions  are 
a  result  of  nonlinear  lift  curve  slope,  specifically  the  fact  that  the  lift  initially  increases  with 
effective  angle  of  attack,  reaches  a  peak  at  45°  and  then  decreases.  The  vehicle  z  velocity 
leads  to  change  in  effective  angle  of  attack  and  thus  one  may  find  multiple  solutions  to  the 
trim  problem. 

Instead  of  prescribing  the  control  variables,  one  can  prescribe  the  trim  state  and  calculate 
the  control  variables.  The  trim  variables  for  hover,  climb  (1  m/s),  and  forward  flight  (1  m/s) 
are  shown  in  Table  5.4.  k0  tilts  the  mean  position  of  the  wings  (majority  of  the  stroke  motion) 
forward  or  backward  of  the  center  of  gravity  of  the  system  providing  pitch  control  (effectively 
the  elevator  when  compared  to  traditional  control  surfaces).  rmaa,  alters  the  magnitude  of 
the  vertical  force  produced  by  the  wings  as  shown  in  Figure  5.9  and  effectively  acts  as  the 
engine,  controlling  forward  or  upward  velocity  depending  on  pitch. 
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Variable 

Trim  Value 

u 

0 

w 

0 

Q 

0 

0 

2.7121e-15 

Kq 

6.9195e-15 

7 max 

1.4784 

(a) 


Variable 

Trim  Value 

u 

-1.7576e-16 

w 

-1 

q 

0 

0 

-1.7576e-16 

K0 

-1.2832e-16 

Tmax 

1.2522 

(b) 


Variable 

Trim  Value 

u 

1.0000 

w 

0.001200 

q 

0 

0 

0.001200 

K  0 

0.003644 

Tmax 

1.5015 

Table  5.4:  Trim  variables  for  prescribed  states:  (a)  Hover  (b)  Climb  (c)  Forward  Flight 


In  order  to  maintain  a  hover  state  it  is  necessary  for  the  vehicle  to  alter  the  lift  magnitude 
from  the  conditions  seen  in  Table  5.3  so  that  it  just  counteracts  the  force  of  gravity  acting 
on  the  body  and  wings.  As  the  hinge  point  is  directly  at  the  center  of  gravity  no  k0  change 
is  necessary  to  balance  the  pitching  moment  and  only  rmax  is  changed,  being  increased  to 
reduce  the  effective  angle  of  attack  and  thereby  the  lift  force  as  seen  in  Table  5.4  (a). 

In  Table  5.4  (b)  the  climb  condition  shows  similar  trim  variables  to  those  of  the  hover 
condition  but  with  a  reduced  Tmax  in  order  to  increase  the  upward  force  as  shown  in  Figure 
5.9. 

The  forward  flight  case  in  (c)  shows  a  u  not  quite  equal  to  one  and  a  w  not  quite  equal 
to  zero  as  those  are  the  body  velocities  and  the  inertial  velocities  were  the  ones  prescribed 
as  described  above.  The  body  is  also  tipped  so  that  the  forcing  vector  can  be  realigned. 
In  this  particular  solution  the  body  is  tipped  back  and  the  forward  and  downward  velocity 
combine  to  create  an  inertial  forward  only  velocity.  In  order  to  allow  the  body  to  move  in 
this  manner  k0  is  used  as  a  pitch  control  to  pitch  the  body  and  Tmax  is  increased  slightly, 
reducing  the  forcing  magnitude  to  the  needed  levels.  While  physically  it  seems  odd  to  tip 
backward  to  move  forward,  the  inertial  forward  velocity  changes  the  angle  of  attack  of  the 
wing  such  that  the  aerodynamic  forcing  still  cancels  the  gravitational  force  appropriately. 
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The  addition  of  the  forward  velocity  also  creates  a  pitch  up  moment  (Figure  5.6)  which 
can  only  be  counteracted  by  the  vehicle  pitching  back  and  using  gravity  to  counteract  the 
moment.  The  tip  back  would  not  be  possible  with  viscous  effects  which  would  force  the 
vehicle  to  tip  forward  to  counteract  the  viscous  drag  (as  discussed  later). 

5.3.3  Stroke- Averaged  Linearized  Stability  About  Trim 


Eigenvalues 

Eigenvalues 

-0.009238  +  0.005345i 

-0.009901 

-0.009238  -  0.005345i 

0.006250  +  0.0091581 

-0.009118 

0.006250  -  0.009158i 

0.006210 

0.01126 

(a) 

(b) 

Eigenvalues 

Eigenvalues 

-0.003819  +  0.008713i 

-0.02662 

-0.003819  -  0.008713i 

-0.009801 

0.01159 

0.001548  +  0.008993i 

0.01174 

0.001548  -  0.0089931 

(d) 


Table  5.5:  Multiple  sets  of  eigenvalues  for  prescribed  kinematics 


Table  5.5  shows  the  eigenvalues  of  the  system  at  the  four  trim  conditions  presented  in  Ta¬ 
ble  5.3.  All  cases  show  instability  about  the  trim  condition  as  indicated  by  the  positive 
eigenvalues. 
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Eigenvalues 

Eigenvalues 

Eigenvalues 

-0.01056 

-0.01190 

-0.01117 

-0.004422 

-0.004867 

-0.001363  +  0.003563i 

0.002210  +  0.003718i 

0.001388  +  0.003279i 

-0.001363  -  0.003563i 

0.002210  -  0.003718i 

0.001388  -  0.003279i 

0.003129 

(b) 


Table  5.6:  Eigenvalues  for  prescribed  states:  (a)  Hover  (b)  Climb  (c)  Forward  Flight 


Table  5.6  show  the  eigenvalues  for  the  hover,  climb  and  forward  flight  trim  cases  presented 
in  Table  5.4.  Again  all  of  these  states  show  instability  about  trim  which  is  supported  by  the 
previous  work  of  Bierling  [1]  and  Richter  [51].  The  transition  between  hover  and  climb  can 
be  examined  through  the  root  locus  plot  shown  in  Figure  5.10  and  the  transition  between 
hover  and  forward  flight  similarly  in  Figure  5.11. 


Figure  5.10:  Root  locus  plot  showing  the  transition  from  hover  to  climb 
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Figure  5.11:  Root  locus  plot  showing  the  transition  from  hover  to  forward  flight 

When  looking  at  the  transition  from  hover  to  climb  all  eigenvalues  maintain  their  sign 
meaning  that  those  that  are  stable  remain  stable  and  those  that  are  unstable  remain  unstable. 
In  the  case  of  hover  transitioning  to  forward  flight,  the  two  oscillatory  eigenvalues  that 
are  unstable  in  hover  move  to  stable  values  for  forward  flight  while  one  of  the  stable  real 
eigenvalues  becomes  unstable. 

5.3.4  Periodic  Trim 

The  actual  periodic  time-varying  system  has  a  periodic  trim.  This  trim  is  calculated  by 
periodic  shooting  as  described  in  Section  4.4.  In  the  case  of  prescribed  kinematics  (kq  = 
0,  Umax  —  0.95944,  r0  =  7t/2,  rmax  =  7t/4,  and  no  deviation)  the  trim  variables  consist  only 
of  the  four  symmetric  states  and  the  resulting  trim  conditions  are  shown  below  in  Figure  5.7. 
The  tables  show  the  trim  variables  which  are  the  initial  values  for  the  states.  The  tables  also 
show  the  stroke-averaged  values  at  trim. 
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State 

Trim  Value 

Average 

u 

0.1037 

2.0063e-08 

w 

-3.4746 

-3.4976 

q 

1.7521 

1.2979e-06 

© 

0.01998 

3.2286e-09 

(a) 


State 

Trim  Value 

Average 

u 

0.3757 

4.8539e-09 

w 

4.5879 

4.5965 

q 

7.2558 

2.8861e-06 

0 

3.1057 

3.1416 

State 

Trim  Value 

Average 

u 

-0.01347 

3.3569e-10 

w 

1.6621 

1.6553 

q 

-0.6694 

-1.7268e-08 

0 

-0.002547 

-1.3051e-10 

(b) 


State 

Trim  Value 

Average 

u 

-0.3542 

6.6106e-09 

w 

-7.09407 

-7.1227 

q 

-1.9779 

1.3349e-07 

0 

3.1394 

3.1416 

(d) 


Table  5.7:  Multiple  sets  of  trim  variables  for  prescribed  kinematics 


Providing  the  initial  conditions  in  Table  5.7  results  in  a  periodic  response  with  average 
values  as  shown.  When  considering  that  the  maximum  states  are  on  the  order  of  10°  these 
average  trim  states  compare  favorably  to  those  found  for  the  stroke-averaged  system  shown 
in  Table  5.3. 

Prescribing  states  and  allowing  two  of  the  kinematic  variables  to  vary  in  order  to  achieve 
hover,  climb,  and  forward  flight  provide  the  results  shown  in  Table  5.8.  Averaging  the  states 
over  one  period  when  providing  these  trim  variables  provides  a  better  comparison  to  the 
stroke-averaged  system  and  results  are  shown  beside  their  initial  value.  As  in  the  prescribed 
kinematic  case  the  averaged  values  (and  the  corresponding  trim  kinematic  variables)  compare 
favorably  with  those  of  the  stroke  averaged  system  shown  in  Table  5.4. 
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Variable 

Trim  Value 

Average 

u 

0.0004158 

-1.2196e-12 

w 

0.01501 

-6.8662e-06 

q 

0.4592 

3.0978e-09 

0 

0.01387 

2.6342e-12 

ft  0 

2.5644e-13 

7 max 

1.4667 

Variable 

Trim  Value 

Average 

u 

0.03585 

7.6309e-10 

w 

-0.9795 

-1.0000 

q 

1.2189 

4.0828e-07 

0 

0.01687 

9.0399e-10 

ft  0 

-1.4761e-12 

7 ~max 

1.2382 

(b) 


Variable 

Trim  Value 

Average 

u 

1.002468 

0.9996 

w 

0.07567 

0.0056 

q 

1.3743 

5.7435e-007 

0 

0.02025 

0.0057 

fto 

0.01006 

^ ~max 

1.4914 

Table  5.8:  Trim  variables  for  prescribed  states:  (a)  Hover  (b)  Climb  (c)  Forward  Flight 


While  the  forward  and  vertical  velocities  listed  in  Table  5.8  (c)  are  not  one  and  zero 
respectively  this  is  because  they  are  body  direction  velocities.  When  they  are  rotated  to  the 
inertial  directions  the  horizontal  velocity  and  vertical  velocity  are  iXs  =  1.0000  and  jig  = 
—4.0172  •  10”4  respectively. 
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5.3.5  Stability  of  the  Time- Varying  System 


Eigenvalues 

Eigenvalues 

-0.008955 

-0.007239 

-0.007726 

0.008161  +  0.02242i 

0.003259  +  0.0061591 

0.008161  -  0.02242i 

0.003259  -  0.0061591 

0.009764 

(b) 


Eigenvalues 

Eigenvalues 

0.001406  +  0.02622i 

-0.009036 

0.001406  -  0.02622i 

-0.0081775  +  0.04753i 

0.003368 

-0.0081775  -  0.04753i 

0.01096 

-0.002184 

(d) 


Table  5.9:  Multiple  sets  of  eigenvalues  for  prescribed  kinematics 


All  of  the  trim  states  in  Table  5.7  are  unstable  as  demonstrated  by  their  positive  eigenvalues 
shown  in  Table  5.9  except  for  the  final  trim  state  which  has  only  negative  eigenvalues.  Com¬ 
paring  these  eigenvalues  to  the  analogous  cases  found  for  the  stroke-averaged  system  (shown 
in  Table  5.5)  it  is  clear  that  the  stroke-averaged  and  time-varying  (periodic)  systems  are  not 
completely  interchangeable.  While  they  do  have  some  similarities  (such  as  trim  conditions 
and  averaged  response)  their  stabilities  are  not  necessarily  the  same.  In  both  cases  they  show 
the  majority  of  trim  conditions  to  be  unstable,  the  only  difference  being  that  the  periodic 
system  shows  one  trim  condition  to  be  stable. 
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Eigenvalues 

Eigenvalues 

Eigenvalues 

-0.008955 

-0.01008 

-0.01541 

-0.007726 

-0.004870 

-0.001014 

0.003259  +  0.006156i 

0.0001577  +  0.002010i 

0.002956  +  0.003488i 

0.003259  -  0.0061561 

0.0001577  -  0.002010i 

0.002956  -  0.003488i 

(b) 


Table  5.10:  Eigenvalues  for  prescribed  states:  (a)  Hover  (b)  Climb  (c)  Forward  Flight 


As  in  the  prescribed  kinematic  case  the  eigenvalues  of  the  time- varying  system  are  on  the 
same  order  as  those  of  the  stroke-averaged  system  but  are  not  similar  otherwise.  This  points 
to  the  fact  that  the  stroke-averaged  system  is  not  a  direct  replacement  for  the  time-varying 
system  despite  the  fact  that  the  flapping  frequency  is  much  higher  than  the  frequency  of  the 
body  dynamics. 


Figure  5.12:  Root  locus  plot  showing  the  transition  from  hover  to  climb 
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Figure  5.13:  Root  locus  plot  showing  the  transition  from  hover  to  forward  flight 


The  transition  from  hover  to  climb  and  hover  to  forward  flight  in  Figures  5.12  and  5.13 
both  show  that  the  system  remains  unstable  through  the  entirety  of  both  transitions.  Also 
those  eigenvalues  that  are  stable  remain  stable  and  those  that  are  unstable  remain  unstable. 
While  the  exact  paths  and  values  are  not  replicated  in  the  stroke  averaged  case  some  of 
the  general  trends  do  appear  in  both  cases.  For  instance,  in  the  case  of  the  transition 
between  hover  and  climb,  the  eigenvalues  with  real  and  imaginary  components  decrease  their 
magnitude  (for  both  the  real  and  imaginary  parts)  in  both  the  stroke-averaged  and  time- 
varying  cases. 

Because  these  trim  conditions  are  unstable  perturbation  from  these  conditions  will  cause 
the  system  to  move  away  from  the  periodic  trim. 

As  an  example  Figure  5.14  shows  the  hover  trim  for  the  nominal  periodic  system  perturbed 
by  one  percent  of  each  value  shown  in  Table  5.8  (a).  The  body  velocities  at  first  stay  near 
the  trim  condition  but  as  the  instability  due  to  the  perturbation  grows  they  quickly  move 
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away  from  the  trim  condition. 


Figure  5.14:  Perturbed 
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5.4  Inflow 

The  first  model  added  for  comparison  is  a  prescribed  inflow  model  as  described  in  Subsection 
3.2.1.  This  constant  inflow  is  based  on  momentum  disk  theory  and  provides  a  constant 
downwash  velocity  (in  the  body  frame).  This  provides  a  first  level  approximation  of  the 
unsteady  portions  of  the  flow  that  would  contribute  to  the  aerodynamics. 
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5.4.1  Periodic  Trim 


Variable 

Trim  Value 

Average 

u 

0.02154 

-4.5330e-ll 

w 

0.01509 

-2.8670e-05 

q 

0.9363 

1.0289e-07 

0 

0.01520 

3.4166e-10 

K0 

7.1690e-13 

T~max 

1.1881 

Variable 

Trim  Value 

Average 

u 

0.04655 

1.9149e-09 

w 

-0.9818 

-1.0000 

q 

1.3176 

5.7236e-07 

0 

0.01674 

2.1190e-09 

ft  0 

6.4085e-ll 

7 max 

0.9849 

(b) 


Variable 

Trim  Value 

Average 

u 

1.01260 

0.9862 

w 

-0.07876 

-0.1627 

q 

2.1133 

1.6284e-06 

0 

-0.1470 

-0.1634 

fto 

-0.01033 

^ ~max 

1.1966 

Table  5.11:  Trim  variables  for  prescribed  states:  (a)  Hover  (b)  Climb  (c)  Forward  Flight 


In  comparing  the  trim  values  for  the  nominal  case  with  a  constant  inflow  model  to  those  of 
the  nominal  case,  some  interesting  similarities  arise.  The  trends  for  all  three  states  are  similar 
in  that  forward  flight  had  the  highest  rotation  magnitude  in  both  cases  while  climb  had  the 
lowest  (other  states  remained  similar  as  well)  which  makes  sense  given  that  the  same  states 
were  prescribed  for  the  both  cases.  It  is  also  interesting  to  note  that  all  of  the  prescribed 
states  showed  lower  values  of  max  rotation  indicating  that  a  geometric  angle  of  attack  was 
needed,  given  the  same  states,  to  counteract  the  decrease  in  angle  of  attack  due  to  prescribed 
inflow.  The  forward  flight  case  here  shows  a  stroke-averaged  tip  forward  rather  than  the  tip 
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backward  of  the  nominal  case  which  is  easier  to  grasp  conceptually. 


5.4.2  Stability  of  the  Time- Varying  System 


Eigenvalues 

Eigenvalues 

Eigenvalues 

-0.01009 

-0.009184 

-0.01767 

-0.003805  +  0.003115i 

-0.006871  +  0.0 1053i 

-0.003195  +  0.006549i 

-0.003805  -  0.003115i 

-0.006871  -  0.010531 

-0.003195  -  0.0065491 

0.002389 

0.005937 

0.008247 

(b) 


Table  5.12:  Eigenvalues  for  prescribed  states:  (a)  Hover  (b)  Climb  (c)  Forward  Flight 


Table  5.12  presents  the  eigenvalues  of  the  system  with  prescribed  inflow.  In  both  Figures 
5.15  and  5.16  those  eigenvalues  that  are  stable  remain  stable  and  those  that  are  unstable 
remain  unstable. 
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Figure  5.15:  Root  locus  plot  showing  the  transition  from  hover  to  climb 


Figure  5.16:  Root  locus  plot  showing  the  transition  from  hover  to  forward  flight 
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Examining  Figure  5.17  shows  that  the  values  for  the  eigenvalues  of  the  various  climbs 
with  inflow  are  near  those  of  climb  without  inflow  at  the  inflow  speed  (-1.23022591896135 
m/s).  This  makes  sense  as  inflow  is  analogous  to  climb  in  many  ways.  The  values  are  not 
going  to  be  exact,  however,  as  the  inflow  is  modeled  to  be  in  the  body  2  direction  at  all  times 
where  as  climb  is  in  the  inertial  z  direction. 


Figure  5.17:  Root  locus  plot  showing  the  transition  from  hover  to  climb  with  a  2  m/s  climb  (w/o 
inflow) 


5.5  Viscous  Drag  Effects 

This  section  adds  a  viscous  drag  model  to  the  nominal  case  as  described  in  Equations  3.33, 
3.34,  and  3.35.  The  assumed  drag  constants  are  Cd(0)  =  0.07  and  Cd( 7t/2)  =  3.06  which 
match  the  hawkmoth  empirical  parameters  provided  by  Berman  [26]. 
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5.5.1  Periodic  Trim 


Variable 

Trim  Value 

Average 

u 

0.02353 

-6.3107e-ll 

w 

0.01507 

-1.5980e-05 

q 

0.9983 

1.8235e-07 

0 

0.01390 

4.1235e-10 

K  0 

2.8213e-ll 

7 max 

1.4673 

Variable 

Trim  Value 

Average 

u 

0.09253 

-1.2566e-08 

w 

-0.9794 

-1.0000 

q 

2.3854 

-5.7279e-07 

0 

0.01340 

2.1154e-09 

K  0 

-8.1659e-10 

7 max 

1.1523 

(b) 


Variable 

Trim  Value 

Average 

u 

1.02438 

0.9893 

w 

-0.06986 

-0.1431 

q 

2.03599 

3.6671e-07 

0 

-0.1286 

-0.1436 

«  0 

-0.04407 

^ ~max 

1.4592 

Table  5.13:  Trim  variables  for  prescribed  states:  (a)  Hover  (b)  Climb  (c)  Forward  Flight 


The  only  significant  change  shown  between  the  nominal  case  trim  variables  and  those  of  the 
nominal  case  with  viscous  drag  are  the  initial  pitch  velocity  values  which  are  significantly 
higher  for  the  nominal  case  with  viscous  drag.  It  is  interesting  to  note  that  even  with 
viscous  drag  the  maximum  rotation  angles  do  not  change  a  great  deal  in  magnitude  (only 
a  few  percent  of  their  nominal  case  values)  to  overcome  the  viscous  effects.  This  indicates 
that  using  the  current  models  assumed  in  this  study,  inflow  (unsteady  effects)  seem  to  have 
a  much  larger  effect  on  the  overall  motion  of  the  body  for  these  three  prescribed  states. 
The  forward  flight  case  exhibits  a  similar  tip  forward  to  the  inflow  case,  which  in  the  case  of 
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viscous  effects  is  necessary.  The  forward  flight  of  the  vehicle  causes  extra  drag  on  the  forward 
stroke  as  compared  to  back  stroke  and  therefore  the  lift  vector  must  be  tipped  forward  to 
counteract  this  force. 

5.5.2  Stability  of  the  Time- Varying  System 


Eigenvalues 

Eigenvalues 

Eigenvalues 

-0.01069 

-0.01357 

-0.01918 

-0.006397  +  0.0081841 

-0.009158  +  0.0153U 

-0.004635  +  0.007753i 

-0.006397  -  0.0081841 

-0.009158  -  0.0153U 

-0.004635  -  0.007753i 

0.007932 

0.007203 

0.01132 

(b) 


Table  5.14:  Eigenvalues  for  prescribed  states:  (a)  Hover  (b)  Climb  (c)  Forward  Flight 


Unlike  in  the  case  of  trim,  viscosity  changes  the  stability  of  the  nominal  model  significantly. 
The  eigenvalues  for  the  nominal  case  with  viscous  drag  show  more  similarity  with  those  of 
the  nominal  case  with  inflow  (Table  5.12)  than  those  of  just  the  nominal  case  (Table  5.10). 
If  the  root  locus  plots  (Figures  5.18  and  5.19)  are  examined  however,  neither  case  shows  the 
same  trends  as  the  nominal  case  with  viscous  drag.  Also  similar  to  the  nominal  case  with 
inflow  (Figures  5.15  and  5.16),  the  nominal  case  with  viscous  drag  shows  a  larger  imaginary 
magnitude  (on  those  eigenvalues  with  imaginary  parts)  than  the  nominal  case  (Figures  5.12 
and  5.13)  alone.  In  the  case  of  moving  from  hover  to  climb  this  magnitude  grows  while 
moving  from  hover  to  forward  flight  sees  a  slight  reduction.  In  both  transitions  from  hover 
those  eigenvalues  that  are  stable  remain  stable  and  those  that  are  unstable  remain  unstable. 
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Figure  5.18:  Root  locus  plot  showing  the  transition  from  hover  to  climb 


Figure  5.19:  Root  locus  plot  showing  the  transition  from  hover  to  forward  flight 


79 


CHAPTER  5.  RESULTS 


5.6  Deviation 


This  section  adds  deviation  from  the  stroke  plane  to  the  kinematic  motion.  The  added 
kinematic  values  are  /i0  =  0,  Umax  =  —0.1744,  and  =  tt/2. 


5.6.1  Periodic  Trim 


Variable 

Trim  Value 

Average 

u 

0.01635 

-8.8938e-ll 

w 

0.03643 

-2.2226e-04 

Q 

1.7920 

3.2974e-07 

0 

0.005028 

1.8819e-10 

K0 

1.7747e-10 

7 ~max 

1.08971 

Variable 

Trim  Value 

Average 

u 

0.03606 

4.9307e-10 

w 

-0.9586 

-1.0000 

Q 

2.2421 

5.3350e-07 

0 

0.01036 

7.1291e-10 

Kq 

4.2675e-ll 

7 max 

0.9557 

(b) 


Variable 

Trim  Value 

Average 

u 

1.01450 

0.9988 

w 

0.05214 

-0.0427 

Q 

2.5756 

2.0669e-06 

0 

-0.03650 

-0.0425 

K0 

0.04563 

^ max 

1.1253 

Table  5.15:  Trim  variables  for  prescribed  states:  (a)  Hover  (b)  Climb  (c)  Forward  Flight 


The  addition  of  deviation  does  make  a  few  changes  to  the  trim  conditions  produced  for 
hover,  climb,  and  forward  flight.  As  seen  in  Table  5.15  the  forward  flight  case  shows  a  tip 
forward  similar  to  that  seen  in  the  cases  with  inflow  and  viscous  effects  added.  Also  the 
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max  rotation  angle  is  decreased  (increasing  the  geometric  angle  of  attack).  The  addition  of 
deviation  removes  body  lift  over  the  majority  of  the  stroke  as  seen  in  Figure  5.20. 


Figure  5.20:  Effect  of  deviation  on  body  forces 


5.6.2  Stability  of  the  Time- Varying  System 


Eigenvalues 

Eigenvalues 

Eigenvalues 

-0.01046 

-0.01038 

-0.01160 

-0.004668 

-0.005489 

-0.003793 

0.006717  +  0.009378i 

0.005612  +  0.005475i 

0.006515  +  0.007638i 

0.006717  -  0.009378i 

0.005612  -  0.005475i 

0.006515  -  0.007638i 

(b) 


Table  5.16:  Eigenvalues  for  prescribed  states:  (a)  Hover  (b)  Climb  (c)  Forward  Flight 

Comparing  the  eigenvalues  for  the  case  with  deviation  (Table  5.16)  to  those  of  the  nominal 
case  (Table  5.10)  shows  the  most  similarity  of  any  of  the  three  effects  so  far  explored.  The 
oscillatory  modes  are  unstable  for  all  three  conditions  in  both  cases  while  the  static  modes 
are  stable.  This  is  clearly  visible  in  the  transition  from  hover  to  climb  and  hover  to  forward 
flight  in  Figures  5.21  and  5.22.  Comparing  these  to  Figures  5.12  and  5.13  shows  similar 
though  not  identical  magnitudes  and  trends. 
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Figure  5.21:  Root  locus  plot  showing  the  transition  from  hover  to  climb 
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Figure  5.22:  Root  locus  plot  showing  the  transition  from  hover  to  forward  flight 
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5.7  Movement  of  the  Hinge  with  Respect  to  the  CG 

5.7.1  Varying  Hinge  Location  in  the  Body  z  Direction 

This  section  alters  the  system  parameters  provided  in  Section  5.1  by  varying  the  hinge 
location  in  the  body  z  (vertical)  direction  from  0  (CG)  to  40  mm  above  the  CG. 

5. 7. 1.1  Periodic  Trim 


Variable 

Trim  Value 

Average 

u 

0.0005014 

7.7938e-ll 

w 

0.01374 

7.8008e-06 

q 

0.2497 

-3.4385e-09 

0 

0.04137 

-8.7338e-10 

«  0 

1.2000e-12 

^ ~max 

1.4301 

Variable 

Trim  Value 

Average 

u 

0.0004158 

-1.2196e-12 

w 

0.01501 

-6.8662e-06 

q 

0.4592 

3.0978e-09 

0 

0.01387 

2.6342e-12 

«  0 

2.50  1  le- 13 

1~max 

1.4667 

(a)  (b) 


Variable 

Trim  Value 

Average 

u 

0.0003320 

-2.0541e-10 

w 

0.01543 

-2.5079e-05 

q 

0.4565 

2.2161e-09 

0 

-0.005092 

-3.1827e-10 

K0 

-2.6800e-14 

7 ~max 

1.45759 

Table  5.17:  Trim  variables  for  prescribed  states  while  varying  the  hinge  location  in  the  body  z 
direction:  (a)  0.04  m  (b)  -0.01  m  (nominal)  (c)  -0.04 
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Examining  Table  5.17  shows  that  moving  the  hinge  in  the  body  z  direction  has  very  little 
effect  on  the  trim  of  the  system,  only  slightly  changing  the  values  of  any  of  the  trim  variables. 

5.7. 1.2  Stability  of  the  Time- Varying  System 


Eigenvalues 

Eigenvalues 

Eigenvalues 

-0.008175 

-0.008955 

-0.008565 

-0.007534 

-0.007726 

-0.007534 

0.001993  +  0.005386i 

0.003259  +  0.0061561 

0.004450  +  0.006626i 

0.001993  -  0.005386i 

0.003259  -  0.006156i 

0.004450  -  0.0066261 

(b) 


Table  5.18:  Eigenvalues  for  prescribed  states  while  varying  the  hinge  location  in  the  body  z  direction: 
(a)  0.04  m  (b)  -0.01  m  (nominal)  (c)  -0.04 


Examining  Table  5.18  and  Figure  5.23  shows  very  little  change  in  the  eigenvalues  as  the 
hinge  vector  changes  in  the  body  z  direction.  One  of  the  static  modes  (the  largest  negative 
magnitude)  increases  in  magnitude  until  the  hinge  vector  is  about  -12  mm  (the  other  static 
mode  begins  moving  to  the  right  before  this  point).  After  that  point  all  of  the  modes  move 
to  the  right. 
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Figure  5.23:  Root  locus  plot  for  varying  the  hinge  location  in  the  body  z  direction 


5.7.2  Varying  Hinge  Location  in  the  Body  x  Direction 

This  section  alters  the  system  parameters  provided  in  Section  5.1  by  varying  the  hinge 
location  in  the  body  x  (horizontal)  direction  from  -  4  mm  to  4  mm  above  the  CG  while 
holding  the  body  z  location  of  the  hinge  constant  at  -10  mm. 
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5. 7. 2.1  Periodic  Trim 


Variable 

Trim  Value 

Average 

u 

0.0004324 

1.1464e-05 

w 

0.01382 

-6.8759e-06 

q 

0.4205 

-2.1598e-07 

© 

0.05374 

0.0400 

«  0 

0.1403 

1 ~max 

1.4671 

Variable 

Trim  Value 

Average 

u 

0.0004158 

-1.2196e-12 

w 

0.01501 

-6.8662e-06 

q 

0.4592 

3.0978e-09 

0 

0.01387 

2.6342e-12 

«  0 

2.5644e-13 

1 ~max 

1.4667 

(b) 


Variable 

Trim  Value 

Average 

u 

0.0003537 

-1.1464e-005 

w 

0.01623 

-6.8755e-006 

q 

0.4878 

-8.1143e-010 

0 

-0.02231 

-0.0400 

K0 

-0.1403 

7 ~max 

1.4671 

Table  5.19:  Trim  variables  for  prescribed  states  while  varying  the  hinge  location  in  the  body  x 
direction:  (a)  -0.004  m  (b)  0  m  (nominal)  (c)  0.004 


Examining  Table  5.19  shows  that  moving  the  hinge  in  the  body  z  direction  has  very  little 
effect  on  the  trim  the  only  large  effect  being  the  shift  of  the  average  pitch  angle  and  stroke 
offset  to  offset  the  shift  in  the  hinge  vector. 
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5. 7. 2. 2  Stability  of  the  Time- Varying  System 


Eigenvalues 

Eigenvalues 

Eigenvalues 

-0.009157 

-0.008955 

-0.009157 

-0.007457 

-0.007726 

-0.007457 

0.003159  +  0.005989i 

0.003259  +  0.006156i 

0.003159  +  0.005989i 

0.003159  -  0.005989i 

0.003259  -  0.006156i 

0.003159  -  0.005989i 

(b) 


Table  5.20:  Eigenvalues  for  prescribed  states  while  varying  the  hinge  location  in  the  body  x  direction: 
(a)  -0.004  m  (b)  0  m  (nominal)  (c)  0.004 


Examining  Table  5.20  and  Figure  5.24  shows  almost  no  change  in  the  eigenvalues  as  the  hinge 
vector  changes  in  the  body  x  direction.  The  largest  changes  are  in  the  static  modes  which 
move  towards  each  other  as  the  hinge  moves  towards  the  CG  from  behind  and  then  back 
to  their  original  values  as  the  hinge  moves  forward  of  the  CG.  Overall  this  indicates  that 
the  movement  of  the  hinge  has  little  effect  on  the  stability  of  the  system  in  hover  though  it 
would  have  a  larger  effect  in  forward  flight  modes  as  this  would  give  the  aerodynamic  forcing 
a  larger  moment  arm  to  act  upon  the  body. 
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Figure  5.24:  Root  locus  plot  for  varying  the  hinge  location  in  the  body  x  direction 
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Conclusion  and  Future  Work 


In  this  thesis  the  influence  of  various  models  on  the  flight  dynamics  of  a  flapping  wing  MAV 
was  investigated.  In  support  of  this  an  aerodynamic  model  and  flight  dynamic  model  were 
developed  in  MATLAB  and  combined  into  a  single  simulation  model. 

The  kinetic  energy  of  the  system  was  derived  using  the  quasi-coordinate  form  of  La¬ 
grange’s  equations  with  rigid  wings  and  fuselage.  The  body  was  assumed  to  be  an  ellipsoid 
while  the  wings  were  rectangular  and  flapped  through  prescribed  kinematics  defined  by  three 
sinusoidal  Euler  angles. 

The  aerodynamic  model  was  based  on  the  Peters  airloads  model  and  allowed  for  the  in¬ 
clusion  of  unsteady  affects  through  the  addition  of  an  inflow  model  and  viscous  drag  through 
an  empirical  model.  These  models  can  be  easily  changed  to  allow  for  more  complicated  wake 
or  viscosity  models.  The  aerodynamic  forcing  was  validated  against  the  quasi-steady  aerody¬ 
namic  model  developed  by  Berman  [26],  including  viscosity  as  implemented  by  Stanford  [3]. 
When  provided  with  matching  parameters  the  models  were  shown  to  be  in  good  agreement. 

For  stability  analysis  a  stroke-averaged  model  and  periodic  model  were  used.  The  stroke- 
averaged  model  averaged  the  forces  and  mass  matrix  over  one  stroke  period.  The  stroke- 
averaged  trim  was  calculated  by  the  Newton  Raphson  method  until  zero-averaged  forcing  was 
achieved  for  prescribed  kinematics  or  flight  condition.  This  was  followed  by  the  calculation 
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of  a  linearized  stability  matrix  about  this  trim  and  from  there  stability. 

The  periodic  system  utilized  periodic  shooting  to  find  trim,  in  this  case  calculating  an  end 
state  after  one  period  of  motion  that  was  equal  to  the  initial  condition.  The  state  transition 
matrix  about  this  trim  then  provides  stability  through  Floquet  analysis.  Comparing  the 
stroke-averaged  and  periodic  systems  indicates  that  there  are  differences  that  make  them 
not  completely  interchangeable.  A  constant  inflow  model  and  viscosity  model  were  added  to 
the  periodic  system  so  that  the  stability  could  be  compared  to  the  nominal  case  (not  including 
either).  While  the  system  remained  unstable  with  the  addition  of  the  models,  it  was  shown 
that  both  have  a  significant  effect  on  stability  (though  not  necessarily  trim)  indicating  that 
for  this  model  (and  most  likely  this  flight  regime)  including  both  of  these  models  is  important. 
Moving  the  hinge  with  respect  to  the  center  of  gravity  had  little  affect  on  hover  stability  but 
with  increased  physical  modeling  fidelity  such  as  body  drag  and  introducing  a  stroke  plane 
angle  offset  from  the  body  angles  the  hinge  location  could  have  a  much  larger  effect. 

For  future  research,  inclusion  of  more  accurate  inflow  models  and  viscosity  models  as  well 
as  more  complicated  physical  models  would  give  more  accurate  comparisons  to  real  world 
applications.  The  addition  of  flexibility  to  the  model  as  well  as  a  body  drag  model  would 
also  further  increase  the  applicability  of  the  model  to  a  real  world  simulation.  A  controller 
could  also  be  designed  to  stabilize  the  system. 
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